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Abstract 

We report quantum Monte Carlo calculations of ground and low-lying excited 
states for nuclei with A < 7 using a realistic Hamiltonian containing the 
Argonne vis two-nucleon and Urbana IX three- nucleon potentials. A detailed 
description of the Green's function Monte Carlo algorithm for systems with 
state-dependent potentials is given and a number of tests of its convergence 
and accuracy are performed. We find that the Hamiltonian being used results 
in ground states of both ^Li and ^Li that are stable against breakup into 
subclusters, but somewhat underbound compared to experiment. We also 
have results for ^He, ^He, and their isobaric analogs. The known excitation 
spectra of all these nuclei are reproduced reasonably well and we predict a 
number of excited states in ^He and ^He. We also present spin-polarized one- 
body and several different two-body density distributions. These are the first 
microscopic calculations that directly produce nuclear shell structure from 
realistic interactions that fit AW scattering data. 

PACS numbers: 21.10.-k, 21.45+v, 21.60.Ka 
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I. INTRODUCTION 



A major goal in nuclear physics is to understand how nuclear binding, stability, and 
structure arise from the underlying interactions between individual nucleons. A large amount 
of empirical information about the nucleon-nucleon scattering problem has been accumulated 
over time, resulting in ever more sophisticated NN potential models. However, for many 
years, it has been feasible to calculate exactly only three- and four-body nuclei with realistic 
NN interactions. Thanks to the ongoing advances in computational resources, particularly 
the advent of massively parallel computers, it is now possible to apply sophisticated quantum 
Monte Carlo methods to the study of light p-shell nuclei, which have a far richer spectrum 
to investigate. 

The p-shell has long been a testing ground for shell model methods |jl|]. Shell-model 
studies of y4 = 4 — 7 nuclei have recently progressed to the stage of large-basis, "no-core" 
calculations using G-matrices derived from the latest AW-interaction models [Q]. Alterna- 
tively, there have been extensive studies of these light nuclei with cluster-cluster models, 
using combinations of AW and cN potentials 0. Our goal here is to calculate properties, in 
particular ground-state energies, of the light p-shell nuclei directly from bare NN and NNN 
interactions, without any intermediate effective interaction. 

Recently we reported |^ first results for the ground states of ^He, ^Li, and ^Be, and the 
first two excited states in ^Li, calculated with the Green's function Monte Carlo (GFMC) 
method for a Hamiltonian containing the Argonne f is NN and Urbana IX NNN potentials. 
Here we present improved and extended calculations for the A = 6 nuclei and the first 
detailed studies of ''He, ''Li, ^Be, and ^B, including all the observed low-lying excitations. In 
this work, it is possible to see for the first time the full splendor of the nuclear shell structure 
emanating directly from a microscopic interaction that fits NN scattering data. 

The Argonne vis model is one of a class of new, highly accurate NN potentials 
that fit both pp and np scattering data up to 350 MeV with a x^/datum near 1. This 
necessarily involves the introduction of charge-independence breaking in the strong force; a 
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complete electromagnetic interaction is also included as an integral part of the model. This 
makes the model useful for studying charge dependence and charge-symmetry breaking in 
nuclei. The AW potential is supplemented by a three-nucleon interaction from the Urbana 
series of NNN potentials , including both long-range two-pion exchange and a short-range 
phenomenological component. The Urbana model IX is adjusted to reproduce the binding 
energy of '^H and give a reasonable saturation density in nuclear matter when used with 
Argonne vis- Details of the Hamiltonian are given in Sec. II. 

The first step in our calculation is the construction of suitable trial functions. Varia- 
tional wave functions based on products of correlated operators have been used successfully 
for and ^He, giving binding energies about 2% above exact Faddeev, hyperspherical 
harmonic, or GFMC solutions 0. We generalize this type of trial function for A = 6,7 nu- 
clei, adding p-wave orbitals and using LS coupling to produce all possible {.P]T) quantum 
states. Parameters in the trial functions are adjusted to minimize the energy expectation 
value, evaluated with Monte Carlo integration, subject to the constraint that rms radii are 
close to the experimental values for ^Li and ^Li. Unfortunately, the best trial functions we 
have been able to build do not give p-shell nuclei stable against breakup into subclusters. 
Nevertheless, these trial functions provide a good starting point for the GFMC calculation. 
The variational wave functions and a brief description of the variational Monte Carlo (VMC) 
calculations are given in Sec. III. 

The GFMC method projects out the exact lowest energy state, \E'o, for a given set 
of quantum numbers, from a suitable trial function, using \E'o = lim^_,oo exp[— (iJ — 
Eq)t\^t- The method has been used with great success in a variety of condensed matter 
problems, and in s-shell nuclei with realistic interactions . Our first calculations for p-shell 
nuclei IQ were made with a short time approximation for the propagation in imaginary time, 
carried out to r = 0.06 MeV~^ and extrapolated to r = oo. In the present work we have 
improved the algorithm by adopting an exact two-body propagator which allows bigger time 
steps, saving significantly on the computational cost. We also have started our calculations 
with better trial functions, which allows more reliable results to be obtained from GFMC 
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propagations that are limited to small r. A detailed discussion of the method as applied to 
reahstic nuclear forces is given in Sec. IV. 

The VMC and GFMC calculations for p-shell nuclei are very computer intensive, and 
would not have been possible without the recent advances in computational power due to the 
advent of massively parallel machines. Section V describes the implementation of the GFMC 
algorithm in a parallel environment, including issues of communication between processors 
and load-balancing. We have also performed a number of tests of the GFMC method, 
including comparisons to other exact calculations for s-shell nuclei, studies of extrapolation 
in r, and sensitivity to the quality of the input trial function. These tests are described in 
Sec. VI. 

Results of our GFMC calculations are presented in Sec. VII. We have obtained energies 
for five states of unique ( J'^; T) in the A = 6 nuclei, and another five states in A = 7 nuclei, 
not counting isobaric analog states. In general we find the nuclei are slightly underbound 
with the present Hamiltonian, but ^Li (^Li) is stable against breakup into a + d {a + t). 
The low-lying excited states are correctly ordered with reasonable excitation energies. The 
VMC energies are found to lie ~ 3 (4.5) MeV above the GFMC results in ^ = 6 (7) nuclei, 
but with very similar excitation energies. We have also used the VMC wave functions to 
perform small-basis diagonalizations of states with the same quantum numbers but different 
symmetries. These calculations optimize the admixtures of different-symmetry contributions 
to the ground state, and also provide estimates for higher-lying excited states with the same 
( J'^; T) quantum numbers. We verify that these states, not all of which have been observed, 
do indeed lie at moderately higher excitations. The VMC spectra are also discussed in 
Sec. VII. 

One- and two-body density distributions from both VMC and GFMC calculations are 
presented in Sec. VIII. These include the densities of spin up/down nucleons in polarized ^Li 
and ^Li. In general, the GFMC densities are slightly more peaked than the input VMC den- 
sities, but overall they are very similar. Unfortunately, both VMC and GFMC calculations 
are not very sensitive to the very long range properties of the wave functions. Therefore 
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it is not yet possible to accurately calculate the quadrupole moments and asymptotic D/S 
ratios with these methods. Finally, we present our conclusions in Sec. IX. 



II. HAMILTONIAN 

Our Hamiltonian includes a nonrelativistic one-body kinetic energy, the Argonne fis 
two-nucleon potential P and the Urbana IX three- nucleon potential @], 

H = Y.K, + Y.v^,+ E y^^k- (2-1) 

i i<j i<j<k 

The kinetic energy operator has charge-independent (CI) and charge-symmetry-breaking 
(CSB) components, the latter due to the difference in proton and neutron masses, 

= Kf^ + Kf'^ ^ + ^)V.^ - ^(^ - ^)r..V? . (2.2) 

4 TUp nin 4 nip rUn 

The Argonne fig potential can be written as a sum of electromagnetic and one-pion- 
exchange terms and a shorter-range phenomenological part, 

v.,=vr+v^+vj^ . (2.3) 

The electromagnetic terms include one- and two-photon-exchange Coulomb interaction, vac- 
uum polarization, Darwin-Foldy, and magnetic moment terms, with appropriate proton and 
neutron form factors: 

v^pp) = Vciipp) + Vc2 + VvP + Vdf + Vmm{pp) , (2.4) 
v^{np) = Vci{np) + VMM{np) , (2.5) 
v^{nn) = VMM{nn) . (2.6) 

The Vmm contain spin-spin, tensor, and spin-orbit components. Detailed expressions for 
these terms, including the form factors, are given in Ref . . These terms should be included 
in calculations aiming for better than 99% accuracy. For example, the contribution of w'^(np) 
to the binding energy of the deuteron is ~ 0.02 MeV, i.e. 1% of the total. 



The one-pion-exchange part of the potential includes the charge-dependent (CD) terms 
due to the difference in neutral and charged pion masses. It can be written in an operator 
format as 

vl = f[^y l^c^ [X,jn ■ r, + X,,T,,] , (2.7) 
where Tjj = Sr^iT^j — ti ■ tj is the isotensor operator and 

X,, = i (Xj + 2x5) , (2.8) 

= I (^S - ^) ' (2-9) 
X™ = [Y{mr)a, ■ aj + T{mr)Sij] . (2.10) 

Here Y{mr) and T{mr) are the normal Yukawa and tensor functions with a cutoff specified 
in Ref. and X^'*^ are calculated with m = m.^± and m^o. 

The one-pion-exchange and the remaining phenomenological part of the potential can be 
written as a sum of eighteen operators, 

^.^■ + ^5= E Mr,,)oi . (2.11) 

p=l,18 

The first fourteen are charge- independent, 

Of='''' = [1, (a, ■ a,), S,,, (L ■ S), L^, L'{a, ■ a,), (L • S)^] ® [1, (r, ■ r,)] , (2.12) 
and the last four, 

0^=''''' = [1, (ai ■ aj), Sij] ® Ti„ {r,i + r,,-) , (2.13) 

break charge independence. We will refer to the potential from the p = 15 — 17 terms as 
v*-^^ and from the p = 18 term as v'"^^ . We note that in the context of isospin symmetry 
the CI, CSB and CD terms are respectively isoscalar, isovector and isotensor. 

The potential was fit directly to the Nijmegen NN scattering data base containing 



1787 pp and 2514 np data in the range — 350 MeV, with a P^r datum of 1.09. It was 
also fit to the nn scattering length measured in d{'K~ ,'-))nn experiments and the deuteron 
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binding energy. There could in principle be more charge-independence-breaking (CIB) terms 
such as L ■ STij or Sijir^i + r^j) but the scattering data are not sufficiently precise to identify 
them at present. 

The Urbana series of three-nucleon potentials is written as a sum of two-pion-exchange 
and shorter-range phenomenological terms, 

V^k = V^^, + V,%. (2.14) 

The two-pion-exchange term can be expressed simply as 

V^^^= Y: A2AXJ',,XJ,}{n.T„r,-n} + C2Ax:,,XJ,][n.T„T,.n] , (2.15) 

cyclic 

where X[j is constructed with an average pion mass, m,r = ^m^o + |m^±. The anticommu- 
tator and commutator terms are denoted by V^^^ and V^j^, respectively, and for the Urbana 
models C2tt = as in the original Fujita-Miyazawa model |TT|. The shorter-range phe- 

nomenological term is given by 

^5 = E V^T\m^n,)T''{m^r,k) . (2.16) 

cyclic 

The parameters for model IX are = -0.0293 MeV and \Jq = 0.0048 MeV. They have 
been determined by fitting the density of nuclear matter and the binding energy of in 
conjunction with the Argonne vis interaction. These values are only slightly different from 
the model VIII values, = -0.028 MeV and Uq = 0.005 MeV, that were adjusted for 
use with the older Argonne fi4 interaction. In principle, the l^^^ can have other terms ||12|| , 



however we need additional data to obtain their strengths; presumably a part of it is due to 



relativistic effects ||T3|-|T5|. 



Direct GFMC calculations with the full interaction (in particular spin-dependent terms 
which involve the square of the momentum operator) have very large statistical errors, for 
reasons that will be discussed in section IV. Also the CIB terms in H are fairly weak and 
therefore can be treated conveniently as a ffist order perturbation. Further, using a wave 
function of good isospin significantly reduces the cost of the calculations. Hence we construct 
the GFMC propagator with a simpler isoscalar Hamiltonian, 



^' = E^r + E<+ E yi^k. (2.17) 

i i<3 i<j<k 



where v',, is defined as 
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E v',{r,j)0!, + v'c{r,j). (2.18) 

p=l,8 



The interaction v^j has only eight operator terms, with operators [1, (cti • ctj), S'jj, (L • S)] iS> 
[1, (ti • Tj)], chosen such that it equals the isoscalar part of the full interaction in all S and P 
waves as well as in the ^Di wave and its coupling to the ^5*1. The strong interaction terms 
are related to the full Vij by 

v[^Vi + f + fvio + fvii + lvi2 + 1^13 + 1^14 , 



v' 



V2 + \vc) + fwio + |f 11 - f f 12 + 1^13 + , 
^^3 + Ivg + f-Uio + fvii - lvi2 + lvi3 + ^Vu , 



^5 = - ^^^13 - ^^^14 , 



Vj^Vt- ^Vg + 1^10 - 5^11 + 1^12 - 1^13 + ^-Vu , 



v' 



^V8 + \vg - hio + ^Vii - hi2 + fvi3 - f V14 . (2.19) 



The isoscalar part of Vci{pp) is also included in H'. We derive it by writing the projector 
for a pair of protons in terms of isoscalar, isovector, and isotensor operators: 

i(l + r,i)(l + r,j) = i(l + In ■ Tj + T,i + T,j + iTij) , (2.20) 
v'cirij) = [ac(A, Z, T) + ^r^ • T^]Vci{pp) , (2.21) 

ac{A, Z, T) = ^^[Z{Z -l) + \A- \T{T + 1)] . (2.22) 

The sum over all pairs of ac{A, Z, T) + ■ tj is just the number of pp pairs in the given 
nucleus. 
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The v^j is a little more attractive than Vij; for example, ^He is overbound by ~ 2 MeV 
with v'^y The expectation value of the difference {vij — v'^j) scales like {Vijk), presumably 
because three-body and higher-order clusters give important contributions to it. Note that 
in and ^He the two-body cluster gives zero contribution to {vij — v'^j) since they are 
identical in low partial waves. We compensate for this tendency towards overbinding by 
using a V-^^ in which the repulsive Uq term of Eq.( p.l6|) has been increased by ~ 30% in the 
H'. This ensures {H') ^ {H) so that the GFMC propagation does not produce excessively 
large densities due to overbinding. The contribution of {H — H') is calculated perturbatively. 



III. VARIATIONAL MONTE CARLO 

The variational method can be used to obtain approximate solutions to the many-body 
Schrodinger equation, = E"^, for a wide range of nuclear systems, including few-body 
nuclei, light closed shell nuclei, nuclear matter, and neutron stars [|T6|. A suitably param- 
eterized wave function, is used to calculate an upper bound to the exact ground-state 
energy, 



Ey = 





H\^v) 







The parameters in are varied to minimize Ey, and the lowest value is taken as the 
approximate ground-state energy. 

Upper bounds to excited states are also obtainable, either from standard VMC calcu- 
lations if they have different quantum numbers from the ground state, or from small-basis 
diagonalizations if they have the same quantum numbers. The corresponding can then 
be used to calculate other properties, such as particle density or electromagnetic moments, 
or it can be used as the starting point for a Green's function Monte Carlo calculation. In 
this section we first describe our ansatz for for the light p-shell nuclei and then briefly 
review how the expectation value, Eq.( p.l| ), is evaluated and the parameters of \l/y are fixed. 
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A. Wave Function 



Our best variational wave function for the nuclei studied here has the form 



1+ E iu^,k+ur,n+T.u, 



LS 



i<j<k 



(3.2) 



where the pair wave function, \l/p, is given by 



i<j 



(3.3) 



The Uij, U^j^ , Uijk, and are noncommuting two- and three-nucleon correlation opera- 

tors, and the 5 is a symmetrization operator. The form of the totally antisymmetric Jastrow 
wave function, \E'j, depends on the nuclear state under investigation. For the s-shell nuclei 
we use the simple form 



II fijk 
i<j<k 



Kj 



\^a{JMTT,)) . 



(3.4) 



Here fci^ij) and f^^j^ are central two- and three-body correlation functions and 

m\\\\))=A\p]p[n]) , 
|$4(0000)) =A\p]pin]ni). 



(3.5) 
(3.6) 



The two-body correlation operator Uij is a sum of spin, isospin, and tensor terms: 



p=2,6 



n fijki^ik,1^jk) 



(3.7) 



while the two-body spin-orbit correlation operator is given by 



(3.8) 



p=7,i 



with 



Of- ^'^ = [1, cTi ■ aj, Sij, (L • S)ij] (g) [1, n ■ Tj 



(3.9) 



The Uij and U^j^ correlations are induced by the two-nucleon interaction. The equations used 
to generate the functions fci^ij) and Up{rij) are given in Ref. [|l7l; they contain a number of 
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variational parameters to be determined by minimizing the energy. The shape parameters 
hsted in Table | seem to have negligible A dependence. Their values are determined by 
minimizing the '^H energy, and are then used for all subsequent A > A calculations. There 
are also a number of parameters that describe the long-range behavior of the correlation 
functions which do vary with A, as well as with the Hamiltonian. Our best values for these 
parameters are given in Table |T|. 

The fiji^, f^ji^, and Uijk are three-nucleon correlations induced by Vy. The first two have 
an operator-independent form: 

fijk = 1 + Qii^ij ■ rifc)(rji ■ rjfc)(rfci ■ Ykj)exp{-qlRijk) , (3.10) 

fijk = 1 - - fifc ■ rjk)exp{-q^Rijk) , (3.11) 

where Rijk = + rik + rjk- The Uijk involve additional spin-isospin operators and are 
somewhat more complicated; they are discussed in Ref. [0]. The U™^ are three-body 
correlations induced by the three-nucleon interaction, which we take in the form suggested 
by perturbation theory: 

t^Sr = E e.^.Jfc(r.„ r,., h^) , (3.12) 

X 

with r = yr, y a scaling parameter, and a strength parameter. Here x = A, C, and R 
refers to the different parts of the NNN potential. With the present Hamiltonian we use the 
three-body parameters given in Table |ITl| for all nuclei studied here. 

The relative importance and cost of computing different elements of the full variational 
wave function are shown in Table |^ for the case of ^Li. The pair wave function \I'p is 
the least expensive to compute, but gives a rather poor energy. The successive addition of 
, U^j^ , and Uijk correlations to make up the full lowers the energy by ~ 2.7 MeV, 
but requires ~ 2.7 times more computation than \E'p. However, about 60% of the energy 
gain can be obtained for only a 20% increase in computation by using the trial function 
defined by 
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i<j<k 



(3.13) 



Here U™' is a truncated TNI correlation from which the commutator term, ec^^fc; has 
been omitted: 



ttTNI _ ~ irA I ^ irR 



(3.14) 



The strength of the anticommutator term is increased to compensate, with ~ l.Se^. From 



Table |I^ we see that this simphfication gets 90% of the gain of adding the full Uf,^\, at 



the cost. The computer time is reduced so significantly because {Xj^.,XJ^} can be written 
as a generalized tensor operator involving the spins of only nucleons i and fc; thus the time 
scales as the number of pairs rather than the number of triples. As discussed below, is 
in fact the most economical starting point for the GFMC calculations. 

The Jastrow wave function for A = Q nuclei is more complicated, as two nucleons must 
be placed in the j9-shell. We use LS coupling to obtain the desired JM value of a given state, 
as suggested in shell-model studies of p-shell nuclei Different possible LS combinations 
lead to multiple components in the Jastrow wave function. We also allow for the possibility 
that the central correlations fdjij) and f^^f^ could depend upon the shells (s or p) occupied 
by the particles and on the LS coupling. The Jastrow wave function is taken as 



nesss TT fssp pssp TT p 
Jijk 11 Jlm5Jlm6 11 Ji 
i<j<k<i l<m<A: 



spp 
n56 



n<4 



n fssir^J) n fsp{rk5)fsp{rk6) E {f3LsfppHr56)\MLSJMTn) 1234:56/ 

j<j<4 



. (3.15) 



fc<4 LS 

The operator A indicates an antisymmetric sum over all possible partitions of the six par- 
ticles into 4 s-shell and 2 p-shell ones. For the two-body correlations we use /ss(r) = fc{r) 
from the ^He wave function, while 



fsp{r) = [asp + bspW{r)]fc{r) + Csp{l - exp[-{r/dsp)^]) , 
fL'ir) = [a,, + b,,Wir)]f,ir) + cji^(l - exp[- (r / d,,)']) . 



(3.16) 
(3.17) 



12 



Here we have supplemented the fdf) with a long-range tail and allowed for a short-range 
modification with a Woods-Saxon factor, W{r) = {1 + exp[(r — i?j)/aj]}^\ The a^p, bsp, 
etc., are variational parameters, whose values are given in Table |V[ For the three-body 
correlations, our best present trial function has f^j^ = f^j^. as in Eg. ( |3.10| ), and f^jj^ = 

Jijk — ^■ 

The LS components of the single-particle wave function are given by: 

|$6(L5JMTT3)i234:56) = | $4(0000) i2340f (i?a5)0f(i?a6) 

A<\h)u,{\t'^)]TT,) . (3.18) 

The 4>p^{Rak) are p-wave solutions of a particle of reduced mass ^miy in an effective a-N 
potential. They are functions of the distance between the center of mass of the a core 
(which contains particles 1-4 in this partition) and nucleon k, and again may be different 
for different LS components. We use a Woods-Saxon potential well: 

VaN{r) = V^'^'ll + expC-^)]-' , (3.19) 

Clp 

where V^^ , Rp, and Op are variational parameters and we allow the depth of the well to vary 
with the LS composition. The values of these parameters are also given in Table 0. The 
wave function is translationally invariant, hence there is no spurious center of mass motion. 
The experimental spectra for A = 6 nuclei |18[ are shown in Fig. |l[ The ground state of 



^He is strong stable, but decays by (3^ emission with a mean life of 807 ms. The first excited 
state is above the threshold for decay to a + n + n and has a width of ~ 1 10 keV; we treat 
it here as a stable state with zero width. In the shell model, the (J'^;T) = (0+; 1) ground 
state of ^He is predominantly a "^^^^Lln] ="'^S[2] state, where we use spectroscopic notation 
to denote the total L and 5* of the state and the Young pattern [n] to indicate the spatial 
symmetry. The (2^^; 1) first excited state is predominantly a ^D[2] state. We allow for a 
possible ^P[ll] admixture in both states, using amplitudes /3oo and Pu in the ground state, 
and (320 and jSu in the excited state. After other parameters in the trial function have been 
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optimized, we make a series of calculations in which the jSis may be different in the left- 
and right-hand-side wave functions to obtain the diagonal and off-diagonal matrix elements 
of the Hamiltonian and the corresponding normalizations and overlaps. We diagonalize the 
resulting 2x2 energy matrices to find the (3ls eigenvectors. The shell model wave functions 
are orthonormal, but the correlated are not. Hence the diagonalizations use generalized 
eigenvalue routines including overlap matrices. We also calculate the position of the three 
complementary 3p[ll] states, with (J^;T) = (2+; 1), 1), and (0+; 1); only one of these 
has been tentatively identified experimentally ||19|. The normalized jSis for these different 



states are given in Table VI 



The spectrum for ^Li contains a (l''';0) ground state that is predominantly ^S[2] in 
character and a triplet of ^D[2] excited states with (3+; 0), (2+; 0), and (1^; 0) components. 
In addition, there are (0"'"; 1) and (2+; 1) excited states that are the isobaric analogs of the 
^He ground and excited states. The ^Li ground state is stable, while the (3^; 0) first excited 
state is above the threshold for breakup into a + d, but is narrow with a width of 24 keV. 
The (0+; 1) second excited state is even narrower, being unable to decay to a + d without 
isospin violation, and thus has a width of only 8 eV. The higher excitations have widths of 
0.5 to 1.7 MeV, but we will treat them here as well-defined states. For the (1"^; 0) ground and 
excited states we mix /?oi, /?2i, and /3io amplitudes by diagonalizing a 3 x 3 energy matrix. 
The latter amplitude corresponds to an unobserved ^P[ll] state, which we also obtain in 
this diagonalization. However, only the (52i amplitude contributes to the (3''"; 0) and (2+; 0) 
excited states. Again, Table ^ gives a summary of the [3ls amplitudes. 

The ground and first excited states in ^Be have the same character as those in ^He, 
except that the Coulomb interaction makes the ground state a resonance, with a width of 
92 keV. Again, we neglect the resonance character of these states in constructing the trial 
function. Most of the numerical results calculated here for the CIB terms of the Hamiltonian 
are obtained by interchanging neutrons and protons in the wave function. This does not 
allow for the changes expected in ^Be compared to ^He due to the Coulomb force, so we 
have also made some calculations adding a Coulomb term V^(r), folded over nuclear form 
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factors, to the VaN^r) used to generate the single-particle 4>p{Rak) functions: 

=2 1 . ... 4 



yi^(r) = Z-{1- exp{-x^)[2 + x^ + — - y-'] 

T Z, 1 y 



21-2 



]exp{-Xp)[2 + Xp + - y^]-^ } 



(3.20) 



Here Xr 



12r/ra, Xp = \J\2rjrp, and y = ra/rp, with the charge radii = 1.65 fm and 



Tp = 0.81 fm. This additional potential term can be used with strength Z = 0, 1, or 2 for 
^He, ^Li, or ^Be, respectively, corresponding to the average Coulomb interaction between 
the a core and a p-shell nucleon. 



The full A = 6 wave function is constructed by acting on the Eg. ( p. 151) , with the 

same Uij, U^^ , Uijk, and U™^ correlations used in '^He. The one exception is that the 
optimal strength of the Uijk correlations is reduced slightly as A increases. In principle, the 
Uij could be generalized to be different according to whether particles i and j are in the s- or 
p-shell, but this would require a larger sum over the different partitions and would increase 
the computational cost by an order of magnitude. 

The Jastrow wave function for A = 7 nuclei is a straightforward extension of Eq. ( p.l5| ) , 
with the added specification of the spatial symmetry [n] of the angular momentum coupling 
of three p-shell nucleons: 



nfSSS TT rssp nssp nssp TT nspp nspp rspp r 
Jijk 11 J ImbJ ImQJ Iml 11 J n56 J n57 J n67 Jl 
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The single-particle wave function is then 
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(3.21) 
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There is an implicit complementary symmetry [n'] for the spin-isospin part of the wave 
function, to preserve the overall antisymmetry, which we do not show here explicitly; for a 
detailed discussion of the symmetry considerations see Appendix IC of Ref. ||20|. The full 
wave function is again built up using Eq.( |3.2| ) with the added definition of the central p-shell 
three-body correlation ff^^ = f^^^. 



The experimental spectra for A = 7 nuclei |]18[ are shown in Fig. In the shell model, 



the lowest states for the T = | nuclei ''Li and ^Be have a predominantly ^P[3] character, split 



2 

into ( J'^; T) = (| , |) ground and (| , |) first excited states. These states are all strong 



2 



stable, the first excited states having mean lives ~ 100 — 200 fs, while the ground state of ''Be 
decays weakly with a mean life of 53 d. Each ground state can also have contributions with 
a mixed spatial symmetry, including ^P[21], ■*P[21], "^Dpl], and ^D[21] components, while 
the first excited state has admixtures of ^P[21], ^P[21], "^Dpl], and ^S[lll] amplitudes. We 
have diagonalized 5x5 matrices for these states. Higher in the spectrum is a predominantly 
F[3] state that splits into (| , ^) and (| , |) pieces. The lower state is mixed with a ^D[21] 
component, while the upper state has ^P[21], ^D[21], and ^D[21] contributions. Again we 
have performed diagonalizations in the PLS[n] amplitudes to project out the lowest (| , |) 
and (| , |) states. 

The diagonalizations confirm that the ground and first excited states are almost pure 
^P[3] and the second and third excited states are almost pure ^F[3]. We have also calculated 
the next excited state of each J as given by our projections to confirm that they lie above 
these first four states. The normalized amplitudes of the lowest two states of each J are 
given in Table [VI I| . These higher excitations include a triplet of predominantly ^P[21] 
states: (| , |), (| , |), and (| , |), and a (| , |) state that is predominantly "^Dpl]. The 
experimental spectrum shows a similar ordering of states, except that the (| , |) comes in 
the midst of the "^Ppl] states, and no second (| , |) state has been identified in this range 
of energy excitation. 

The spectrum for T = ^ states in A = 7 nuclei is also shown in Fig. 0. The (| ; |) 
ground state for ''He is 0.44 MeV above the threshold for breakup into ^He+ra with a width 
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of 160 keV. The isobaric analogs have widths of 260, 320, and 1200 keV for ^Li, ^Be, and ^B, 
respectively. This state can have contributions from ^P[21], ^D[21], and ^S[lll] amplitudes, 
and we again diagonalize a 3 x 3 matrix to evaluate the PlsIu] components, as shown in 
Table |V14 We also calculate the first three excited states: a (| ; |) ^P[21] state, a (| ; |) 
^D[21] state, and the second mixed (| ; |) state. None of these excited states have been 
experimentally identified. 



B. Energy Evaluation 

The energy expectation value of Eg. ( ^31) is evaluated using Monte Carlo integration. 



A detailed technical description of the methods used here can be found in Refs. ||17| , |21| , |22 
Monte Carlo sampling is done both in configuration space and in the order of operators in the 
symmetrized product of Eq. (|3.3| ) by following a Metropolis random walk. The expectation 
value for an operator O is given by 

The subscripts p and q specify the order of operators in the left and right hand side wave 
functions, while the integration runs over the particle coordinates R = (ri, r2, . . . , r^). This 
multidimensional integration is facilitated by introducing a probability distribution, H^pg(R), 
such that 



(0) = ■ (3.23) 





^1{R)0^,{R)/W,,{R) 


; w,,{R) 


J2p,q I dR 


'^l{R)^g{R)/W„{R) 


W,,iR) 



Ep,JrfRK(R)* '^^ 

This probability distribution is taken to be 

W,,{R) = |Re((vl>t,_p(R)vl/p,,(R)))| , (3.25) 

which is constructed from the pair wave function, \I^p, but with only one operator order of 
the symmetrized product. This probability distribution is much less expensive to compute 
than the full wave function of Eq. ( ^.21 ) with its spin-orbit and operator-dependent three-body 
correlations, but it typically has a norm within 1-2% of the full wave function. 
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Expectation values have a statistical error which can be estimated by the standard de- 
viation cr: 



a 



21 1/2 



(3.26) 



N -I 

where N is the number of statistically independent samples. Block averaging schemes can 
be used to estimate the autocorrelation times and determine the statistical error. 

The wave function \l/ can be represented by an array of 2^ x (^) complex numbers, 



vl/(R)=^^„(R)|a) 



(3.27) 



where the '?/'a(R) are the coefficients of each state \a) with specific third components of spin 
and isospin. This gives arrays with 96, 960, 1280, 2688, and 4480 elements for "^He, ^He, 
^Li, ^He, and ^Li, respectively. The spin, isospin, and tensor operators Of^"^'^ contained in 
the two-body correlation operator Uij, and in the Hamiltonian are sparse matrices in this 
basis. For forces that are largely charge-independent, as is the case here, we can replace this 
charge- conserving basis with an isospin-conserving basis that has N[A, T) = 2^ x I{A, T) 
components, where 



I{A,T) 



2T + 1 

:A + T+1 



\ 2 



A 

hA + T 



(3.28) 



This reduces the number of array elements to 32, 576, 320, 1792, and 1792 for the cases 
given above - a significant savings. In practice, the Xj ■ tj operator is more expensive to 
evaluate in this basis, but the overall savings in computation are still large. 

Expectation values of the kinetic energy and spin-orbit potential require the computation 
of first derivatives and diagonal second derivatives of the wave function. These are obtained 
by evaluating the wave function at 6A slightly shifted positions of the coordinates R and 
taking finite differences, as discussed in Ref. |T^. Potential terms quadratic in L require 
mixed second derivatives, which can be obtained by additional wave function evaluations 
and finite differences. A rotation trick can be used to reduce the number of additional 



locations at which the wave function must be evaluated 23 



As a check on the correctness of our Monte Carlo integration, we have evaluated the 
energy expectation value (H) for the deuteron using the exact wave function as input, and 
match the energy to better than 1 keV. We have made similar calculations for the triton 



using a Faddeev wave function, as discussed in Ref. JT^, and obtained agreement with 
independent Faddeev calculations at the 10-20 keV level. For the much more complicated 
A = 6,7 wave functions, we also evaluate the expectation values (J^) and ( J^) to verify that 
they truly have the specified quantum numbers. A third check is made on the antisymmetry 
of the Jastrow wave function by evaluating, at an initial randomized position. 



(3.29) 



where P^j'"'^ are the space, spin, and isospin exchange operators. This value should be 
exactly zero for an antisymmetric wave function, and it is in fact less than 10^^ for each 
pair of particles in each nuclear state that we study. 

A major problem arises in minimizing the variational energy for p-shell nuclei using the 
above wave functions: there is no variational minimum that gives reasonable rms radii. For 
example, the variational energy for ^Li is slightly more bound than for ^He, but is not more 
bound than for separated ^He and nuclei, so the wave function is not stable against 
breakup into a + d subclusters. Consequently, the energy can be lowered toward the sum 
of ^He and energies by making the wave function more and more diffuse. Such a diffuse 
wave function would not be useful for computing other nuclear properties, or as a starting 
point for the GFMC calculation (see Sec. VI below), so we constrain our search for optimal 
variational parameters by requiring the resulting point proton rms radius, r^, to be close to 
the experimental values for ^Li and ''Li ground states. For ^He and ''Be ground states, and 
all the excited or resonant states, there are no experimental measurements of the charge 
radii. To avoid introducing too many additional parameter values, we construct these wave 
functions by making minimal changes to the ^Li and ''Li wave functions, with the added 
requirement that the excited states should not have smaller radii than the ground states. 

For A = 6 nuclei, we begin by selecting parameters to minimize the energy of the ^Li 
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^S[2] component (the dominant part of the ground state) subject to the constraint that 
Tp ~ 2.4 — 2.5 fm. For the other components, only the depth of the single-particle well, 
V^^, of Eq.(pJ9D, and the tail, c^^^, of the p-shell pair correlation function of Eq. (|3.17| ) are 



varied, as shown in Table [V|. The well depth for the ^D[2] states is decreased to get the rms 
radius of the (3~^; 0) excited state larger than the ground state. The tail is increased for the 
mixed-symmetry ^P[ll] state for the same reason. For ^He we use the same parameters as in 
^Li for the corresponding ^S[2], ^D[2], and ^P[ll] states. The only other difference between 
^Li and ^He wave functions is that we may turn on the a-N Coulomb potential of Eq.( p.20[) 



when generating the single-particle radial functions (pp. Finally, the diagonalizations are 
made to determine the I3ls mixing coefficients of Table 0. 

A similar procedure is followed for the A = 7 nuclei. Parameters are selected for the 
dominant ^P[3] state in ^Li subject to the constraint that Tp ~ 2.2 — 2.3 fm. The well 
depth is reduced for the ^F[3] states, and the tail is increased for all the mixed symmetry 
states. Afterwards the (3 is diagonalization is carried out. Since ^Be is a mirror nucleus, it 
has the same wave function as ''Li, aside from changing the a-N Coulomb potential. The 
^He and ground states are isobaric analogs to mixed symmetry states in ^Li, so they use 
corresponding parameters. 

Shell model lore tells us that the lowest state of any given {J'^;T) will be the state with 
maximal spatial symmetry and smallest L that can be formed from the allowed couplings, 
e.g., the ^S[2] ground state in ^Li or the ^P[3] ground state in '^Li. For the purposes of 
obtaining a variational upper bound and a GFMC starting point, we could settle for a 
constructed using only that LS[n\ component. However, by using all the allowed compo- 
nents, we can gain a significant amount of energy in some cases and, as is be discussed below, 
this gain persists in our GFMC propagations. For the A = 6 nuclei, the diagonalizations 
for the (0"'"; 1), (1"'";0), and (2''"; 1) states improve the lower state by 0.25 to 0.5 MeV. In 
the first four ''Li T = | states, the mixing is much less, and improvements are at most 0.15 
MeV. However, for '^He, there is a gain of 0.75 MeV, probably because there are two states 
of identical symmetry that only differ by 1 in L. 
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The diagonalizations have the additional benefit that we can predict where the next 
higher excited state of each ( J'^; T) hes. This allows us to confirm that the Hamiltonian is 
not predicting any unusually low excitations that are not observed experimentally; e.g., the 
second (| ; |) and (| ; |) states in ''Li do not appear below the first (| ; |) and (| ; |) 
states. 



IV. GREEN'S FUNCTION MONTE CARLO 

The aim of the GFMC method is to project out the exact lowest energy state, \&0; asso- 
ciated with a chosen set of quantum numbers, from an approximation to that state. The 
method used here is essentially identical to that used previously to calculate nuclei with A< 
6 0], with the exception that we have now incorporated the exact two-body propagator in 
the imaginary-time propagation. In this section we describe the algorithm in some detail, 
in particular relating it to algorithms commonly used for scalar interactions. For simplicity 
of notation we will not make the distinction between H' and H (and their respective com- 
ponents) that was introduced with Eq.( |2.17D ; the reader will want to remember that we in 



fact use the simpler H' in our GFMC propagator. 

GFMC projects out the lowest energy ground state using = lim^^oo exp[— — 
Eo)t]'$t- The eigenvalue Eo is calculated exactly while other expectation values are generally 
calculated neglecting terms of order |\l/o — ^I/tP and higher. In contrast, the error in the 
variational energy, Ey, is of order |\E'o — ^Tp, and other expectation values calculated with 

have errors of order |\I'o — ^^tI- 

We use the of Eq. ( |3.13|) as our initial trial function and define the propagated wave 
function \1/(t) as 

^(r) = e-(^-^°)"^r; (4.1) 

obviously \E'(r = 0) = and \&(r oo) = \l/o. Introducing a small time step. At, 
T = nAr, gives 
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^-{H-Eo)At 



^T- (4.2) 



The \I/(t) is represented by a vector function of R using Eq.( ^.27D , and the Green's function, 
6*0/3 (R, R') is a matrix function of R and R' in spin-isospin space, defined as 

Gap{K,K') = (R,a|e-(^-^")^"|R',/3). (4.3) 

It is calculated with leading errors of order (Ar)'^ as discussed in Sec. IV.2. The errors in 
the full calculation are determined by the difference between the (artificial) Hamiltonian for 
which the propagator is exact and H. This difference is of order (Ar)^, and Ar is chosen to 
be small enough that this total error is negligible. Omitting spin-isospin indices for brevity, 
\E'(R„, r) is given by 

^(R„, t) = J G(R„, R„_i) ■ ■ ■ G(Ri, Ro)^T(Ro)c^Rn-i ■ ■ ■ Ric^Ro- (4.4) 
The mixed expectation value of an operator O is defined as: 



(O) 



Mixed 





0|*(r)) 




a/(r)) 



/ dP„4(R„)0G'(R„, R„_i) ■ ■ ■ G(Ri, Ro)^t(Ro) 
/ c?P„4(R„)G(R„, R„_i) ■ • • G(Ri, Ro)^t(Ro) 



(4.5) 



where P„ = Rq, Ri, ■ ■ ■ , Rn denotes the 'path', and 

dPn = dRodRi ■ ■ ■ dKn . (4.6) 

In GFMC, the integral over the paths is carried out stochastically. Generally, the required 
expectation values are calculated approximately from the variational and mixed expec- 
tation values. Let 

^(r) = ^r + 5^(r). (4.7) 

Retaining only the terms of order 5\E'(r), we obtain 

{^{t)\0\^{t)) 



(0(r)) 



(nrmr)) ' 

{0{r))M^.ed + [{0{T))M^.ed " {0)t] , (4.8) 
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where 

io)T = • (4-9) 





0\^t) 







More accurate evaluations of (0(r)) are possible, essentially by measuring the observable 
at the mid-point of the path. However, such estimates require a propagation twice as long as 
the mixed estimate. Since we are limited in the present calculations to a total propagation 
time of 0.06MeV~^, we use the approximation ([4.8| ). 

An important exception to the above is the energy, Eq given by {H{t — ^ oo)). The 
{H{t)) Mixed can be re-expressed as (Ref. p5| ) 



{H{T))M^.ed = (^^|e-(H-£;oW2e-(H-i?o)r/2|^^) ' (4-10) 

since the propagator exp[— (if — Eq)t\ commutes with the Hamiltonian. Thus {H{t)) Mixed 
approaches Eq in the limit r oo, and furthermore, being an expectation value of it 
obeys the variational principle 

{H{T))Mixed > Eo . (4.11) 

If a simpler H' is being used to construct the GFMC propagator, then these equations apply 
to {H'ij)), and {{H — H')) must be evaluated using Eq. (^^ ). 

Given these expressions, two basic elements are required for any GFMC calculation. The 
first element is the choice of short-time propagator exp[— (if — Eq)At\ and the second is a 
method for sampling the paths. We discuss each of these elements in turn. 

A. The Short-Time Propagator 

The short-time propagator should allow as large a time step Ar as possible, since the 
total computational time for propagation is proportional to 1/Ar. Earlier calculations 
0J2^J^ used the propagator obtained from the Feynman formulae. Ignoring three-nucleon 
interaction terms in H, it is given by 
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-HAt 



-KAt 



i<j 



+ C(Ar3) . 



(4.12) 



Note that it is useful to symmetrize the product of 6"""'^^^/^ when [vij,Vjk\ 7^ 0, in order 
to reduce the error per iteration. The nuclear Vij has a repulsive core of order GeV. The 
main error in the above propagator comes from terms in e~^^'^ having multiple Vij, like 
VijTvij{AT)^ for example, which can become large when particles i and j are very close. 
In order to make them negligible a rather small Ar ~ 0.1 GeV~^ is used with the above 
propagator. The matrix elements of the propagator are given by: 



5 J^e-^.j(r»j)AT/2 



i<j 



Kj 



-v,,(r' )At/2 



(4.13) 



G'o(R,R') = (R|e-^^"|R') 



m 



27rrAr 



3A 



exp 



fR-R 



/\2 



2n^AT/m 



(4.14) 



However, it is well known from the studies bulk helium atoms that including the 
exact two-body propagator allows much larger time steps. This short-time propagator is 



(jq,^(R, R') — Go(-R') -R-') ('^1 
where Qij is the exact two-body propagator, 



I] J 



_ -H,jAT\ I 



(4.15) 



(4.16) 



m 



and g^^ij is the free two-body propagator. 



go,ij[^ij, r 



27rrAr 



exp 



(4.17) 



(4.18) 



where /i = m/2 is the reduced mass. All terms containing any number of the same Vij and 
K are treated exactly in this propagator, as we have included the imaginary-time equivalent 
of the full two-body scattering amplitude. It still has errors of order (At)^, however they 
are from commutators of terms like VijTvikiAr)^ which become large only when both pairs 
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ij and ik are close. Since this is a rare occurrence, a five times larger time step At ~ 0.5 
GeV^^ can be used for the present studies of light nuclei. In the case of bound states of 
helium atoms a ~ 30 times larger time step can be used with the propagator ( [4.15| ) than 
with ( |4.12| ) presumably because the inter-atomic potentials have a relatively harder core, 
and they commute with each other. 

Finally, including the three-body forces and the Eo in Eq. (^4.3|) , the complete propagator 
is given by 



G«,(R, R') = e^°^-Go(R, R') exp[- + ^.f.(R'))^] 



(a|/3(R)|7)(7l 



\6){5\hiR')\/3) 



(4.19) 



At 



(4.20) 



The exponential of Vj~1 is expanded to first order in At thus, there are additional error terms 
of the form V^lV^^J,f.,{ATy . However, they have negligible effect since V^l has a magnitude 
of only a few MeV. It was verified that the results for ^He do not show any change, outside 
of statistical errors, when At is decreased from 0.5 GeV~^. 



B. Calculation of gij 

The pair propagator g^j is a matrix in the two-body spin-isospin space, and obeys the 
equation 



[-^ + H^,] (7.,(r,r';r)=0. 



(4.21) 



As the Hamiltonian naturally decomposes into eigenstates of the two-body spin and isospin, 
so does the propagator gij. In addition, it obeys the convolution equation 



g,,ir, r'; t + t') = J dh"g,,{v, r"; r)^?,,(r", r'; t') 



(4.22) 



with the initial condition 
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g,,{r,r';T = 0) = S{r,r') . (4.23) 

To calculate gij, we use the techniques developed by Schmidt and Lee [^] for scalar 
interactions. The basic idea is to use the convolution equation, Eq. (^4.22|) , to write Qij as a 
product over N steps: 

N 

QijiTN, ro; Ar) = \{ g^j{r,, ri_i; e), (4.24) 

i=l 

with e = At/N and an implied integration over intermediate points. If we use a symmetric 
expression for the short-time propagator gij{e) such as 

g,,{r, r'; e) = e-^-(r)^/2^o.,(r, r'; e)e-^^^^'^'^^/' , (4.25) 

the errors for gij{AT) contain only even powers of starting with l/N"^. ||28| By evaluating 
gij{r, r'; Ar) for several values of N (and consequently e) and extrapolating to e — >• the gij 
can be calculated with high (~ 10 digit) accuracy. 

Evaluation of gij can be carried out in various ways. We have chosen to expand the prop- 
agator in partial waves denoted by JM, TT^, S, and L, thus replacing the three-dimensional 
integral in Eq.( |4.22| ) with many one-dimensional integrals. The two-nucleon interaction Vij 
has a simple form, Vjj^'g{r), in these partial waves. The interaction is diagonal (L = L') in 
S = and 1 waves with L = J, and it couples the S = 1,L, L' = J±1 waves. 

The gij{r, r'; Ar) is written as a sum over partial waves, 

,.,(r, r'; Ar) = ^ E E XT,T.yr,sir) ''^'^"'"J' y\Ur')xy. , (4-26) 

JM TT^ SLU 

where Xt,t^ denote isospin states, and yfisi^) standard spin- angle functions that depend 
upon spins and the directions r and r'. In uncoupled channels the partial wave propagators 
g^Tsi^i^''i ^'^) scalar functions of the magnitudes r and r', while for coupled channels 
they are 2x2 matrix functions with L, L' = J±l. The partial wave propagators obey the 
one-dimensional convolution equation 




(4.27) 
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where 



gQ .ij is the free propagator in partial waves with angular momentum L, 



9o,ij 



X 



(4.28) 




with 



(4.29) 



X = 



2h\ 



These expressions are used to obtain 5'jt5('^, '?') for several values of e, extrapolating 



these 'channel' propagators very accurately, because they must be summed to reproduce a 
Gaussian falloff in the angular variables that is present in the full propagator. Fast Fourier 
transforms are used to switch between momentum and coordinate space, where the kinetic 
and potential terms, respectively, are diagonal and can be trivially exponentiated. 

We then sum over partial waves to obtain the full two-body propagator. If we were 
to include only the (physical) anti-symmetric two-body channels, the complete two-body 
propagator would also be anti-symmetric, and hence for small r the propagator would have 
two peaks, one near the original point and another (with a minus sign in symmetric spin- 
isospin states), near the point corresponding to the interchange of the particles. In principle 
we could use this propagator by sampling paths with an arbitrary permutation at each 
step, and perhaps cancel some noise arising from unphysical symmetric states. However the 
propagation distance is governed by the Gaussian behavior of Go, which is much shorter- 
ranged than the average pair separation. Hence, any cancellation would be very small. 
Instead, we simply use the Argonne v'^ potential in unphysical states, and include them in 
the propagator. In essence this corresponds to treating the particles as Boltzmann particles 
for purposes of the propagator. Since one always computes overlaps with completely anti- 
symmetric states, this is perfectly acceptable. 

One also has complete freedom to choose an arbitrary interaction in the unphysical 
channels, but the present choice retains the property of a positive definite Green's function 
in spin-singlet channels. The propagator in 5" = states is positive definite, since it is for 



to e 



until a specific error tolerance has been reached. It is important to obtain 
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small e and the convolution form of the Eq. ([4.22|) preserves this property. The choice of 
the Boltzmann propagator allows us to simply sample Gaussians centered on the identity 
permutation when choosing the paths. 

It is also important to include many partial waves in the calculation. The starting 
gij{r,r'; e) is a narrow Gaussian of width [Ah^e/m) which is < 0.1 fm for e < 0.1 GeV~^. 
Hence a large number of partial waves are required to reproduce it accurately. The propa- 
gator in all J < 55 states is calculated from Eq. (|4.27| ); beyond this we use simple approx- 
imations including the analytically known propagators for free particles. Keeping a much 
smaller set of partial waves would yield the same answer in an exact quadrature, however it 
can dramatically increase the statistical error in a Monte Carlo calculation. For example, the 
positive definite property described above is recovered numerically only for large numbers 
of partial waves. 

The terms having and (L- operators are not presently included in the propagator. 
These terms, like others that depend quadratically on the relative momentum between the 
interacting particles, represent changes in the mass of the particles due to interactions. 
They can in principle be included by using appropriate effective masses in the kinetic energy 
propagator gQ^ij [^]. Unfortunately, there is a strong spin-isospin dependence in the and 
{L ■ S)"^ interactions, which then makes the go,ij spin-isospin dependent. Attempts to use 
them have generally led to large statistical errors. For this reason the propagator uses the 
approximate f g interaction operator. 

Calculating this propagator is computationally intensive. Therefore, prior to the GFMC 
calculation, we sum the propagator over partial waves and store the full sum on a grid. 
Storing the gjTsi^", t"'', At) is impractical, both because of memory requirements and the fact 
that summing over waves for each r and r' would be computationally expensive. For a spin- 
independent interaction, the propagator gij would depend only upon the two magnitudes r 
and r' and the angle cos{6) = r • r' between them. Here, though, there is also a dependence 
upon the spin quantization axis. Rotational symmetry allows one to calculate the spin- 
isospin components of gij{r,r') for any r and r' by simple SU3 spin rotations and values 
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of Qij on a grid of initial points r = (0,0,2;) and final points r' = {x',0,z'). The x — z 
plane is chosen because the yjis ^^^^ there. In addition, the fact that the propagator 
is Hermitian allows us to store only the values for z > z'. In the z-direction we take an 
evenly spaced grid of 0.02 fm extending up to 6 fm. Beyond 6 fm the Vij is weak and it 
is sufficient to use Eq.( [4.15| ) to calculate Qij. The propagator falls off approximately as a 
Gaussian, Eq. (^4.18|) , with range parameter 2h^ Ar/m ~ 0.3 fm for Ar = 0.5 GeV^"*^. Thus 
the x' and z — z' grids have maximum values of ~ 0.9 fm and are nonuniform. 



C. Sampling of the Paths 

The remaining task in a GFMC calculation is to sample a set of paths; in order to 
maintain a reasonable statistical error we sum explicitly over all spin-isospin states of the 
system for each path. To choose the paths we follow as closely as possible the standard 
practice for scalar interactions, as we have done in previous work P,p|, |26| , p9| . In this section 
we compare the standard method with that used for nuclear systems. 

The integrals in Eq.([4.5|) for {0{t)) Mixed are carried out stochastically using a relative 
probability function P(P) to sample the paths. Each path consists of a set of n steps, where 
each step contains a sample of 3A particle coordinates, as well as sets of operator orders used 
to sample the symmetrization operators S for the pair operators in the trial wave function, 
Eq.( |3.13D , and the propagator, Eq. (|4.19|) . The ensemble of the sampled paths is denoted by 



{P}; and contains Np paths. For each path P we define 

iVp = 4(RjOG'(R„,R„_i) ■ ■ ■G(Ri,Ro)^T(Ro)/i'(P) , (4.30) 



Dp = ^5.(R„)G(R„,R„_i) • • ■G(Ri,Ro)^t(Ro)/P(P) , (4.31) 



and 



(0(r))M^.ed = (E ^p)/(E ^P) ' (4-32) 
{P} {P} 

with a statistical error determined by the correlated variance of A^p and Dp and proportional 
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to 1/yjNp. The relative probability function -P(P), should be chosen to minimize this 
statistical error. 

Different schemes for sampling the paths with probability -P(P) are possible. In finite- 
temperature simulations, one typically retains the entire history of the path and uses a 
Metropolis scheme to sample them |^|. For zero-temperature simulations, however, it is 



generally more efficient to sample the paths through a branching random walk. Points along 
the path are generated iteratively through an importance-sampling procedure. Only the 
amplitudes of the propagated wave function and the accumulated weight of the path need 
be retained for each configuration. We discuss the algorithm used for nuclear spin-isospin 
dependent interactions after first describing the algorithm for spin-independent interactions. 

For scalar interactions and real Hamiltonians, the particles can be assigned specific spin 
states in \E'7-(R) which never change during propagation, G(R., R') is a real, positive function 
with finite norm, and "^t, and consequently ^E'(t), can be chosen as a real scalar function. 
The -P(P) is commonly taken as 

1 



^(p) = n 



i=l,n 



/(Ri)G'(Ri, R; 



i-1 



/(Ro)|^t(Ro)| . (4.33) 



-^(Rj-l; 

The importance function /(R) is used in sampling and hence should be positive definite, it 
is often taken to be the magnitude of the trial wave function, 

/(R) = |^t(R)| • (4.34) 

The initial configurations are sampled from /(Ro)|\E't(R-o)|- The quantity in brackets in 
Eq.( [4.33|) is referred to as the importance-sampled Green's function Gj, 

1 



G/(Ri, Rj-i^ 



/(Ri)G(Rj, R' 



■i-i. 



(4.35) 



/(R.-i) 

The probability of the path P(P) depends implicitly upon all of the steps in the path, but is 
decomposed into an initial weight /(Ro)|\&r(R-o)|, times a product of weights for each step. 
Using this -P(P) in our expressions for A^p and Dp we arrive at 

/(R„) |^t(Ro)| ^ ' 
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- /(R„) |vI..(Ro)| • ^'-''^ 
In the ideal case of a Bose ground state, \l/r(R) is positive for all R; choosing /(R) = \E't(R) 

yields = 1 with zero variance and the variance of A^p is acceptable for many interesting 

operators. In particular, if \I't'(R) is close to the ground state of H, then 

~ ^° ' ^^-^^^ 
and the A^p for {H{t)) Mixed will have a small variance. Many properties of Bose liquid 
and solid ^He and its drops have been studied with GFMC ||27| , |30|j3T[| using this probability 
density. 

In contrast, the wave functions of simple Fermi systems have domains of positive and 
negative signs separated by nodal surfaces. The importance function I given in Eq.( [4.34|) 



must be increased slightly near the nodal surfaces to allow diffusion between the domains. 
When the path P crosses a nodal surface its Dp and A^p change sign. At small r, few paths 
are long enough to cross nodal surfaces and the variance is small. As r increases, many 
paths cross nodal surfaces, the variance increases and the average value of Dp decreases. 
This problem is called 'the Fermion sign problem', and it limits the maximum value of r up 
to which the state can be propagated. P2| , |5B| Generally the calculations are continued until 
the statistical error increases beyond an 'acceptable' point; we iterate until r = 0.06MeV^^ 
for most nuclei studied here. 

For local spin-isospin independent interactions, ^H, ^He, and ^He nuclei would be com- 
pletely spatially symmetric and no sign problem would exist. For more realistic local interac- 
tions, the dominant spatially-symmetric component of the wave function and the relatively 
large excitation energies imply that the sign problem is not very significant for three- and 



four-body nuclei |]3^. However, it does limit the propagation of states with A > A which 
must have nodal surfaces as required by antisymmetry, and also all calculations with non- 
local interactions. 

Implementing the algorithm to sample the paths is straightforward. Choosing /(R) = 
|\Efr(R)|, the initial (r = 0) configuration Rq for each path is obtained, as in VMC, by 
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sampling ^f.(R) using the Metropolis method. The subsequent configurations Rj, at r = 
zAr, are obtained sequentially from Rj_i, by iterating with the importance-sampled Green's 
function G/, 

/(R,)^(R,) = J G,(R.,R.-i)/(R.-i)^(R^-i)c^Ri-i. (4.39) 

This equation is the importance-sampling generalization of the iterative form of Eq.( [4.4|) , 

^(R,) = J G(R„R,_i)^(R,_i)rfR,_i. (4.40) 

Eq.( |4.39| ) describes the evolution of the density /(Rj)|\l'(Rj)| with r = iAr, hence the 
configurations Rj are distributed with this density. The propagation is entirely in terms of 
distinguishable 'Boltzmann' particles; the Fermi or Bose character of the system is retained 
only at the two ends of the walk through the statistics of the initial and final trial wave 
functions. 

Up to this point, we have assumed that we can sample points along the path directly from 
Gj, but typically this is not possible. One must sample from an approximate Gi(Ri,Ili^i) 
and then use the weighting and branching techniques discussed below to create paths with 
probability proportional to the product of the (^/(Rj, Ri_i). If points are sampled from an 
approximate Gj, it is convenient to define a weight 

mR.R.-i) = £4i'^ (4.41) 

as the ratio of the full Gj to the approximate Gj. Simply choosing paths with P(P) = 
nj=i,n G'/(Ri, Rj--i)/(Ro)|^E'r(Ro)| would modify expressions for the the numerator and de- 
nominator, Eqs.( |4.36| - }4.37| ), by multiplying the contribution of each path by the product of 
the W: 

iVp = iy(P)^^^^^, (4.42) 
^ ' /(RJ |M/t(Ro)| ' ^ ^ 
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where 

W{P)= n W^(Ri,Ri-i) . (4.44) 



4=1, n 



As a trivial example, one could sample the free-particle propagator (G/ = Gq), and the 
weights W would be the ratio of final to initial importance functions times the ratio of 
interacting to free-particle propagators. Such a scheme, however, is woefully inefficient. As 
the path length increases, so do the fluctuations in the W(P), and the branching techniques 
discussed below must eventually be used to control them. 

For an efficient and unbiased calculation, it can be very important to choose a G to 
minimize fluctuations in the weights W introduced at each step. For scalar problems, one 
typically samples a shifted Gaussian, where the shift is related to the logarithmic derivative 
of the trial wave function. This can be used to perform importance sampling accurate to 



second order in Ar, (for a review, see Ref. |^5[|) and hence essentially set Gj = Gj 



In the nuclear case, though, the wave function consists of many spin-isospin amplitudes, 
and a more complex sampling scheme is required. For illustrative purposes, we describe 
the scalar equivalent of our sampling method, although for scalar interactions it is not as 
efficient as sampling a shifted Gaussian. The free propagator, Go(R'jRi-i)5 can be easily 
sampled. A number of points, R^-, j = l,nsamp are obtained by sampling Go(Rj) Rj-i)- 
These points should be chosen in a correlated manner to reduce fluctuations. Anticipating 
requirements for the non-scalar case, we define an approximate scalar importance-sampled 
Green's function Gf(Rj, Rj_i). The primary requirements are that Gf is fast to compute, 
that it is positive, and that it approximates Gj; for the scalar case one could simply choose 

Gf = Gj. 

For each of the risamp points, we calculate Gf(R',Rj_i). The Rj is picked from the set 
with probability proportional to Gf (R^-, Rj_i)/Go(Rj-, Ri-i). This procedure imphcitly 
defines a G/, and requires a weight 

G/(Ri, Rj-i) 



1 ^ G'/(Rj7Rj-ly 

nsarap j=l,„,_p G'o(R;-, Ri_l) 
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Gf(R„R,_i)- 



(4.45) 



Only the variance, and therefore the statistical sampling error in the calculation depends 
on Usamp- When Ugamp = 1; R-j = R-'i a^d the vector Rj — Rj_i can be in any direction 
since Go depends only upon (R'^^ — Rj_i)^. In this case, the weights W = Gj/Gq can differ 
significantly from unity and add to the variance. Indeed the growth estimate of the energy, 
obtained from the difference between unity and the ratio of new to old weights, will have 
an infinite variance in the limit At — 0. In the present calculations, we consider only two 
points H'l and Rg = 2Rj_i — R'^^ symmetric about Ri_i. The leading gradient contribution 
along R'^ — Rj_i, in the expansion: 

is thus cancelled up to order At and the variations in 1^(P) are reduced significantly. 

Nevertheless, the weights of different paths used in computing expectation values, 
Eqs. (|4. 4^4^.431) will eventually diverge. This divergence yields an increasing statistical error. 



as the contribution of only a few paths will dominate the others. Consequently, branching 
techniques are required to control the fluctuations in the relative contributions of different 
paths. In branching, the configurations are redistributed every few time steps by keeping 
n\i unit weight copies of each configuration where 

ni, = znt{W{P) + Cr) , (4.47) 

Cr is a random number between and 1, and int denotes the (truncated) integer part. 
The W^(P) of the resulting configurations are then set to one in order to account for the 
branching process. This branching technique, in effect, forces the paths to be sampled from 
the product of Gj rather than G/. 

On average, the expectation value of any path is reproduced correctly using this tech- 
nique. However, the computation is much more efficient as configurations with small weights 
are more likely be discarded while configurations with large weights are replicated. In this 
way, statistical noise is reduced by keeping an adequate population of contributing configu- 
rations. 

34 



The algorithm used for nuclear GFMC, in which there is a strong spin-isospin dependence 
to the interaction, is a generalization of the procedure described above. Here wave functions 
must be regarded as vectors in spin-isospin space and the (^(R, R') as a matrix, however, 



the relative probability of the paths -P(P) must remain a scalar. Following Refs. p6|,p|,|29 
we define an importance-sampled Green's function Gi as well as an approximate Gj. Just 
as in the scalar case, the approximate Gi is used for sampling points in the path, while the 
ratios Gi / Gj define weights W which are used in branching. 

In order to introduce importance sampling, we first define a scalar function / of the trial 
and GFMC, Eq.(|3), wave functions: 

/[vI/r(R,),vI/,(P,)] = |^4^jR,)vl/,„(P,)| +e^|v^^,jR,)vl>,,(P,)| , (4.48) 

a a 

where \E'r is the trial wave function, \E'j is the i'th iteration of the GFMC wave function, 
which depends implicitly upon the path Pj, and a denotes the spin-isospin components. 

This definition of the importance function differs slightly from the scalar case, which only 
involved the trial function \E't. The first term simply measures the magnitude of the overlap 
of the wave functions, while the second, with a small coefficient e (~ 0.01) ensures a positive 
definite importance function to allow diffusion across nodal surfaces. In this definition of / 
as well as the remaining discussion, we suppress the sampling of the pair orders in the wave 
function and the propagator. 

The importance-sampled Green's function Gj can then be defined as the ratio of the 
importance functions after one iteration of the full propagator G'(Rj,Rj_i), 

rrR R J[vi/HR.),^.(R.)] 

where 

^,(R,) = ^(R,, R,_i)^,_i(R,_i) . (4.50) 



For scalar interactions, this definition of Gj is equivalent to Eq.( ^.35D . Here, the importance 



function is defined from all the amplitudes of the trial and GFMC wave functions, and the 
effects of the propagator are included in the importance function. 
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To perform a calculation, the initial configurations are sampled from /(Po), which is 
defined by inserting \l/j=o = ^^t in Eq.( |4.48| ) above. For speed, the VMC calculations use 
the simple importance function Wpq(R) defined in Eg . [S .251 . Hence we introduce the ratio of 
the two importance functions as an initial weight and perform a branching step immediately. 
This procedure results in a population drawn from /(Pq)- 

Sampling steps from Gj gives us 

= /[*r(Rn),^n(Rn)] , (4.52) 

and hence estimates of observables as the ratio N/ D where 

n.=i,„G(R„R,_i)^o(Ro) 



p 



/[^r(R„),^„(R„)] 



^ _ Vi;^kRn)a=l,nG'(R»R.-l)^o(Ro) 

/[^T(Rn),*n(Rn)] ' ^ ' ^ 



For scalar interactions, setting e = in the definition of /, Eg. ( [4. 481) , we recover Egs. ( [4.3€ - 



OTl) . 

Again, though, we cannot sample from Gi directly. We must sample from a Gj and 
introduce weights and branching to correctly get paths sampled from the products of Gj. 
The procedure is exactly as described previously, although here it is important to introduce 
an approximate Gf . The Gf is a spin-independent function, and hence is much easier to 
compute than Gi, which involves all the spin-isospin states of the system. The present 
algorithm reguires us to compute only a single full propagator per iteration, and the full 
trial wave function only after several iterations. 

The scalar importance function Gf(R',Rj_i) is again used to implicitly define Gj, and 
construct weights M//(Ri, Rj-i). It contains scalar approximations to the dominant physics 
present in the propagator and the trial wave function, 

Gf(R,R') = |vi/,(R)|G^(R,R')^-^ , (4.54) 

where G'^(R, R') is obtained from an approximate spin-isospin independent interaction: 
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vs{rij) = ^[vcCSo,rij) +VcCSi,rij)] (4.55) 
and the Feynman approximation 

G^(R,R') = exp{-Y,vs{r,,)Ar/2)Go{n,n')eM-T.^s{<^)^T/2) . (4.56) 

i<j i<j 

This propagator uses the average of the central potentials in the important S-waves, and, 
like the true G(R, R'), is small at small Vij , preventing the configurations from having 
small inter-particle distances inside the repulsive core range. Similarly, the approximate 
importance sampling in Gj is governed by the function Eq. (|3.4| ), which can be used 
ClS db simple approximation to I^E'tI- 

With these definitions, a step in the propagation is the same as in the scalar case. It 
begins with sampling Ugamp correlated points R^ from the free-particle propagator. Then 
Gj is evaluated for each possible step, and we choose Rj from them. The weight W is then 
computed as the ratio of importance functions divided by the sampling probability G as in 
Eq.(^. 

Since we typically do not compute observables after each step, and fiuctuations in the 
weights are not significant after a single step, it is not necessary to compute the importance 
function / (and hence the trial wave function) at every step. We perform branching after 
every second step, with the weights computed from the product of intermediate W. This 
product involves only the Gj and overlaps of the wave functions at the final and initial 
steps. Hence, for intermediate steps we must compute the full G acting on the GFMC wave 
function, but not necessarily the full trial wave function. Again, after branching, the weight 
of each path is set to W^(P) = 1. 

At this point we can reconstruct the estimates of any observables. After branching, the 
A'p and Dp obtained with this -P(P) are given by Eg. ( [4.531) . Ignoring e, the variance of Dp 



is mostly due to the Fermion sign problem, while that of A^p is tolerable, particularly when 
O = H and is close to the desired eigenfunction of H. For an exact \E'j', we regain the 
exact ground state energy with zero variance. 
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V. COMPUTATIONAL METHODS 



The complicated nature of the nuclear interaction and the computational complexity of 
the calculations presented here require high performance computing. In the past, vector 
supercomputers were used for the first *^He and ^He GFMC calculations. As we have stated 
before, the size of the wave function vector grows exponentially with the number of nucleons 
and the number of matrix operations grows with the number of pairs. In making the step 
from four- to six-body calculations at least an order of magnitude increase in computational 
performance was required. The clear means of achieving this performance goal was parallel 
computation. 

A frequent method of achieving performance gains in Monte Carlo calculations is to 
distribute the configurations over several processors and let each processor carry out its own 
independent Monte Carlo calculation. Such an "embarrassingly parallel" implementation 
is sufficient for simple calculations in which each processor can handle a calculation with 
a minimum acceptable number of configurations in a reasonable amount of time. For the 
seven-body systems considered here this is not the case. 

The heart of the VMC calculation is the Metropolis algorithm which is an inherently 
serial algorithm. Since the bulk of the work in our variational calculations lies in the energy 
expectation value, the straightforward division of labor is to have one master processor 
perform the Metropolis walk while several slave processors calculate the energy and other 
expectation values for the configurations that the master generates. The number of slave 
processors that can be efficiently used is the ratio of the CPU time needed for expectation 
values to that needed to walk from one configuration to the next. We find that typically 50 
processors can be used efficiently in a ^Li VMC calculation. 

Implementing the GFMC algorithm on a parallel architecture provides some special chal- 
lenges. There is no clear division of labor as in the VMC calculation and the number of 
configurations can change throughout the calculation. An embarrassingly parallel implemen- 
tation could work for ^Li, since one processor on current machines is capable of handling 
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enough ^Li configurations (several thousands) for an independent GFMC calculation. For 
^Li, on the other hand, only hundreds of configurations could be propagated on each pro- 
cessor to achieve an acceptable turnaround time for an independent calculation. With such 
a small configuration set on each node, the population fiuctuations on each processor would 
leave some processors with few configurations and others with too many. To avoid such an 
inefficient use of resources, periodic load balancing between processors is required. 

In our implementation, the initial configurations are generated and written to disk in a 
random walk that uses only two processors. Typically 50,000 configurations are generated. 
These are then used in one or more subsequent GFMC calculations. At the start of the 
GFMC calculation, the master processor reads the configurations and distributes them to the 
slave processors. It is then responsible for collecting and averaging energy expectation values, 
and determining load-balancing distributions. Each of the slave processors is responsible for 
a block of configurations. The slaves perform propagation and branching for this block of 
configurations. At selected values of r (typically every 20 steps) they save the configurations 
in local lists for subsequent energy calculations. During load balancing (which typically is 
done every 10 steps) each slave reports its current load to the master, which then instructs 
each over-loaded slave to send its excess configurations to under-loaded slaves. In this way, 
all slaves have, within tolerances, the same number of configurations. The master must 
receive load information from all the slaves before it can determine the redistribution. In 
order that slaves not remain idle while waiting for the redistribution information, they 
compute energies for the configurations stored in their local lists. When a slave completes 
a block of energy calculations for a particular time step, they are sent back to the master. 
After all of the energy results for a time step are received by the master, they are averaged. 
This program structure scales well with the number of processors since there are no major 
communication bottlenecks in the course of the calculation and load balancing keeps the 
slaves somewhat synchronized. Calculations with up to 50 processors on the Argonne IBM 
SP show no degradation in efficiency; typically the slaves are idle less than 5% of the time 
and most of this idle time occurs at the end of the calculation. 
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The current version of our GFMC program is written using Fortran 90 and makes use 
of the MPI message-passing hbrary. On IBM SPl processors, we achieve 40% (33%) of 
the theoretical speed in ^Li (^Li) calculations. On IBM SP2 wide nodes we get 45% for 
both cases. The better efficiency can be ascribed to the larger cache on the wide nodes. 
Table |V111| shows the performance of our GFMC program on the Argonne IBM SP (using 
SPl nodes) and the Cornell IBM SP2 using wide nodes. The table gives the CPU times for a 
single propagation step and a single energy evaluation and the total time required for 50,000 
configurations propagated for 120 time steps (Ar=0.0005 MeV^^, Tfinai=0-06 MeV~^) with 
load balancing every 10 time steps and energy evaluations every 20 steps; this time includes 
slave idle time and average effects of configuration number growth. It is based on results 
using 20-40 slaves. However, this total time does not include the master time. As one can 
see, the time grows by roughly a factor of 10 from ''He to ^Li to ^Li. The total computational 
effort approximately scales as the product of the wave function size N{A, T) = 2^ x I{A, T) 
[Eq. ( p. 281 )] and the number of pairs P = ^A{A — 1). This rapid increase in computation 
time will be a serious obstacle to extending these calculations beyond eight-body nuclei on 
presently available computers. 

Our earlier 6-body calculations |^ were about 30 times slower than those we now make. 
The increased speed is approximately attributable to 1) using the exact two-body propagator 
(factor of 3), 2) dropping the ec^ji term in the (factor 2), 3) improvements in the 
calculation of the wave functions (factor 3), and 4) reduction of processor idle time (factor 
1.4). 

VI. ACCURACY OF GFMC 

In this section we consider several aspects of our calculations that could introduce sys- 
tematic errors in the GFMC results, and attempt to place limits on these errors. There 
are two major sources of error: due to the fermion sign problem, the GFMC propagation 
cannot be extended to arbitrary large imaginary time and thus admixtures of low-energy 
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excitations in the trial wave function will not be fully removed, and the GFMC propagation 
must be done with a different Hamiltonian from the desired one. We also investigate the 
effects of the time step size. 

Figure ^ shows the statistical errors as a function of imaginary time for calculations 
of {H) for various nuclei using 50,000 initial configurations. The errors grow dramatically 
with increasing A due to the increasingly poorer quality of the ^t- In these calculations, 
{'^t\H\^{t)) is evaluated by having H act entirely to the left, therefore if \I/t is an eigen- 
function of if, the sampling errors will be nil as per Eq. (^4.5|) . For the p-shell nuclei, the 
errors increase exponentially with r. This is the well-known fermion sign problem; the s- 



shell nuclei suffer much less from this problem . This exponential error growth places an 
effective limit of rjj„ai=0.06 MeV~^ on our calculations for A = 7 nuclei, which means that 
admixtures of states of excitation energy less than ~12 MeV in '^t-, will be damped by less 
than 50%. 

This led us in Ref. to attempt to extrapolate the computed E{t) to r = oo by fitting 
them with 

where Eq is the extrapolated energy, and the strengths of contaminating states in are 
approximated with a few at excitation energies E*. Figure^ shows such fits made for ^He. 
Because ^He calculations are so inexpensive, we were able to make precise calculations of 
E{t) for many values of r and thus determine the short-time behavior of E{t) using 200,000 
to 740,000 configurations. The solid curve is a fit with Eq = —28.335 MeV, excitation 
energies E* = 20.2, 341, and 1477 MeV, and corresponding = 0.0062, 0.0018, and 
0.00046. The lowest O"*" excitation of ^He is at 20.2 MeV and this energy was not varied 
in the fit. The of the fit is 19 for 31 E{t) (25 degrees of freedom), so the E{t) are 
not statistically independent. We have not attempted to estimate the correlations between 
the energies at different times. The increases by 1 when Eq is changed by +0.02 or 
—0.03. The dashed curve shows a fit without the 20.2 MeV excitation; it gives = 23 and 
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Eq = —28.28 MeV. For most of the other GFMC calculations reported in this paper, we 
did not compute E{t) at the many r < 0.1 MeV^^ used in these fits. Therefore we made 
several fits to the E{t) for r > 0.1. A fit using E^ = 20.2 MeV and one adjustable E* 
gives Eq = —28.33(3) with = 14.6 (11 degrees of freedom), while a fit with just one E* 
results in Eq = — 28.33l;;'j'2, -^i = 30, and = 16.0. Finally, the dash-dot line and dotted 
lines shows the average of the E{t) for 0.04 < r < 0.1 : —28.300(15). It appears that in 
this most favorable case, with high statistics, high first excited state, and large maximum 
r, we can see that including the first excited state improves the extrapolation marginally. 
However, the extrapolated Eq is not significantly lower than a simple average of the E{t) 
for 0.04 < r < 0.1. 

Figure ^ shows the E{t) and fits made for the ground state of ^Li. The values 
for r > 0.06 MeV^^ were computed with 200,000 initial configurations, those for r = 
0, 0.01, 0.06 MeV~^ have 280,000 configurations, while those for the other small r have 
only 50,000 configurations. The energy at very small r is influenced by admixtures of very 
high energy states in These have little effect on the E{t > 0.1 MeV~^, therefore we 
make fits to E{t) only for r > 0.01. The dashed curve is a fit to the E{t) for 0.01 < r < 0.06, 
which is the range that is available for the other p-shell nuclear states in this paper. The 
extrapolated energy is Eq = — 31.56to;5o MeV, where the indicated errors correspond to 
increasing by 1. This fit was made using a single excitation energy, E^ = 36 MeV. The first 
1"*" excited state of ^Li is at 5.65 MeV. A single-energy fit constrained to this energy gives 
large x^- Two-energy fits with one energy constrained to 5.65 MeV have a very flat x^(^o) 
from which useful values of Eq can not be extracted. The solid curve shows a single-energy 
fit made to the E(t) up to 0.1 MeV~^ available for this state ; it gives Eq = — 31.38l!o;i8. 
We see that including data up to 0.1 MeV^^ reduces the error in Eq by about a factor of 
two. Finally the solid line with dotted error bars is the average of the r = 0.04, 0.05, and 
0.06 values, denoted by Eav Its value, -31.25(11) MeV, is formally an upper bound for Eq 
and is above the extrapolated Eq by only one standard deviation. 

Because of the difficulties in making useful extrapolations in r, it is important to under- 
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stand contaminations in particularly from low-excitation-energy states which will not be 
fully filtered out by r = 0.06 MeV~^. We have made several calculations of the ground-state 
of ^Li to study the effects of changes in on the GFMC E{t). Figure ^ shows the effects of 
removing some of the noncentral correlations in "^t] the solid circles are from a calculation 
with the full and are the same as in Fig. ^. The open squares were computed by using 
the simpler \['p of Eq.( p.3|) . This makes the energy at r = worse by ~ 1.7 MeV. However 
by r = 0.01, the GFMC has fully corrected for this defect and thereafter the differences 
are just statistical fluctuations. Hence removing Uf^j^^ from enhances the admixtures 



of excitations > 250 MeV. Calculations without the UT^^^ would be about 20% faster than 



) 

ijk 



full calculations, but the poorer quality of without U^j^ increases the statistical errors 
at large r by about 20%, and thus requires 40% more configurations to get the same error. 
Therefore it is not economical to drop the U™^ from "^t- 

The open diamonds in Fig. ^ come from a much more drastic approximation of Here 
both the and tensor components Uij have been omitted, resulting in a four-operator 

wave function. In such a wave function, the dominant tensor components of the two-body 
potential have zero expectation value and the energy at r = is -|-41 MeV. It is completely 
corrected by r = 0.03; again the rate of correction indicates excitation energies ~ 250 MeV. 
The statistical errors from such a bad are much larger. 

These two tests indicate that defects in the non-central parts of the correlation, which 
have been the subject of much optimization in VMC studies, are easily corrected by the 
GFMC. Deficiencies in the one-body part of \E'r present more of a problem. As is discussed in 
Sec. HI, the for ^He has two symmetry components: [2] and [11]; the optimal amplitudes 
for these (see Table [V^ are 0.967 and —0.253, respectively. The solid circles in Fig. |^ 



show GFMC energies from a using these components. The open squares show results 
computed using a with just the [11] component; the E{t = 0) obtained with such a 
wave function is 4.5 MeV higher than that obtained with the best \E'r; this corresponds 
to the 5 MeV excitation energy of the dominantly [11] state. However, because this error 
is entirely due to a low energy excitation, the GFMC makes very little improvement by 
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r = 0.05 MeV ^. A less radical case is shown by the open diamonds which correspond to a 
with amplitudes of +:^, — for the two states. The E{t) starts out 1.6 MeV above 
that of the best ^r; the GFMC reduces this to ~ 1.3 MeV at r = 0.06 MeV"^ The sohd 
curve is a single-energy fit to these results; the fitted excitation energy is 36 MeV and the 
extrapolated Eq = —26.7 is well above the Eav = —27.64(14) from the best Fits with 
two excitations, one constrained to 5.0 MeV, give an essentially fiat ^ind are not useful 
to extract the Eq. 

A similar situation arose in our first GFMC calculation for the ground state of ''He. This 
was made with just the dominant ^P[21] component. A subsequent calculation using the 
three components listed in Table [VI 1| lowered all the E{t) by ~ 0.8 MeV; the first | excited 



state is at only 2 MeV so we would have to propagate 10 times further for the GFMC to 
substantially correct this error. 

It is important that the \I/t have the correct admixtures of different symmetry, and 
other low-lying states. Otherwise the GFMC results will only be upper bounds to the exact 
eigenenergies. For this reason the used in this work are obtained by diagonalizing the 
Hamiltonian in the small bases of low-energy shell-model states as discussed in Sec. IV. Such 
diagonalizations can be used to ascertain if an improved will infiuence the GFMC results 
obtained with \E't- The Hamiltonian should be diagonalized between "^'j, and \1/t, taking 
into account their nonorthogonality. If the difference in the eigenvalues is large compared 
to 1/Tfinai the results will not be influenced; if it is small, the superposition corresponding 
to the lowest eigenvalue must be used. As examples of this, we made such diagonalizations 
for the cases studied above, in which \l/p or just a four-operator were used. These 
gave excitation energies of ~ 700 MeV and ~ 300 MeV, respectively, which are in good 
agreement with the observed E{t). Our best \I'y contains U^j^ and Uijk correlations omitted 
from the due to computational costs. Diagonalizing the Hamiltonian within and 
shows that these correlations admix states with excitation energies of ~ 1000 and 300 MeV, 
respectively. Thus they can be safely left out of the 

As is discussed in Sec. HI, the VMC calculations for p-shell nuclei do not have a local 
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variational minimum for reasonable rms radii. Therefore the variational searches were con- 
strained to have radii close to the experimental values, if such values are known. To study 
the sensitivity of the GFMC results to this assumption, we have made a number of GFMC 
calculations of the ^Li ground state using that have different rms radii. These \1/t were 
made by changing the depth (Vp) and radius (Rp) of the Woods-Saxon well used to make the 
p-wave orbitals (see Eq. ( p.l9| )); thus the ^He core was not directly modified. Figure || shows 
the evolution of the rms radius with r for these calculations; the solid circles correspond to 
the \1/t used in the rest of this paper. We see that the GFMC basically makes no change 
to the radii, even though they span a range of almost 30% (there may be some indication 
that the smallest radii are increased at small r). This is probably because completely sepa- 
rating the deuteron from the ^He core corresponds to only a 1.5 MeV excitation. Figure |^ 
shows the GFMC energies from these calculations as a function of the GFMC radii (both 
are averages of the 0.04 < r < 0.06 MeV~^ values), and the corresponding expectation 
values. The variational energies obtained with decrease monotonically with increasing 
rms radius, but the GFMC energies show a weak minimum; the very large radii yield higher 
GFMC energies and thus can be variationally ruled out. The curve is a parabolic fit to the 
five GFMC energies with smallest rms radii; the minimum is at 2.44 fm. However the curve 
is very flat and the uncertainty in the location of the minimum is at least 0.1 fm. Thus even 
when 50,000 to 280,000 initial configurations are used for each point, it is difficult to extract 
the ^Li rms radius. Given this, it is not possible at present to reliably study the radii of 
A = 7 nuclei. 

The GFMC propagator used in these calculations has two possible sources of error. The 
first is that the time step, Ar = 0.0005 MeV^^, might be too large. We have checked this by 
two calculations. For we made calculations to r = 0.06 MeV~^ using both Ar = 0.0005 
and 0.00025 MeV'^ these were different by 0.016(14) MeV. For ^He we made calculations to 
r = 0.01 MeV"^ using both Ar = 0.0005 and 0.0001 MeV'^ these were different by 0.02(6) 
MeV. The statistical errors in our GFMC energies for p-shell nuclei are all > 0.1 MeV; thus 
these time-step effects are neghgible. 
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The second possible source is that we cannot use a propagator for the full Hamiltonian, 
if, that we are interested in. Rather we must use the H' of Eq. ( |2.17| ) and compute {H — H') 
perturbatively. Kamada and Glockle have estimated for that evaluating (vij) in wave 
functions appropriate for v'^j underestimates the binding energy by < 20 keV; scaling this 
by the total two-body potential energy gives < 50 keV in ^He and < 90 keV in ^Li. We 



constructed a He that is optimized for v'^j + Vijk and used it to compute this gives 
60(15) keV less binding than our best ^He variational wave function in agreement with the 
above estimates. 

Using a propagator for an H' that gives more binding than H can introduce small errors 
in the determination of densities and radii. The more tightly bound eigenstate of H' is 
likely to have a smaller radius. The radii of '^He can be calculated more accurately, and 
we have studied their sensitivity to various propagators. The rms radius for optimized 
for H is 1.482(3) fm. A GFMC calculation using an H' with no v'(j{rij) and V-jj. = Vijk 
(see Eq. (|2T7| )) gives {H - H') = 2.40(3) MeV and an rms radius of 1.418(4). However 
using the ^^(^ii) and I.SUq in V-jf^ results in {H — H') = .03(2) MeV and an rms radius of 
1.446(3) fm. Presumably the later value is more correct, while the former is too small due 
to the overbinding. 

As is discussed in Sec. IV, the GFMC directly computes mixed estimates {0{t)) Mixed- 
Except for H' and operators that commute with it, these must be corrected to obtain the 
desired (0(r)); we use Eq. (^I8| ) to achieve this. Consequently, the expectation values of 
the individual energy components, such as fC'"^, v[p and V^j^^, which have errors of order 
|vl>g — \l>Tp, do not sum to the correct total energy. Indeed, there must be a collective error 
in these individual terms equal to the total difference between the GFMC {H') Mixed and 
the VMC {H')t- This is illustrated in Table ff^ for the case of ®Li, where the difference 



{H')Mixed — {H')t is —4.4 MeV, and the sum of the individual (O) is an additional —4.4 
MeV lower than {H'). In this case, the individual corrections are comparable in magnitude 
to the collective error, but small compared to the total expectation values. 

Aside from our own work in Ref. [Q, there are no published calculations of p-shell nuclei 
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using realistic interactions such as those used here. However we can compare to previous 
values for the s-shell nuclei. There are accurate Faddeev and projected hyperspherical 
harmonics (PHH) calculations of '^H for the Argonne f 14 with no Vijk- For this Hamiltonian 
we find a GFMC energy of —7.670(8) MeV which is in good agreement with the previous 
results of -7.670 (Faddeev/R 0), -7.680 (Faddeev/Q and -7.683 (PHH A 

PHH result for Argonne vis with Urbana IX has recently been computed [^: -8.475. It is in 



good agreement with our values of -8.455(8) obtained with Ar=0.5 GeV ^ and —8.471(12) 
obtained with Ar=0.25 GeV"^ 

There are also several other calculations of ^He with Argonne f 14 without Vijk- Figure 
p!0| shows the GFMC E{t) for this case. Because there is no Vijk, "we multiplied the v'^j in H' 
by 0.994 and included the isoscalar so that {H — H') ~ 0. The average, shown by the line 
and dotted error range, of the last few E{t) is —24.227(31) MeV. A calculation using the full 
v'ij, without isoscalar v'^, in H', gives —24.230(31), even though in this case {H — H') = 1.5 



MeV. These results are in excellent agreement with the older GFMC calculation of Ref. |^ 
(up-pointing triangle) which was made with a completely independent program that uses 
the short-time propagator of Eq. ( |4.12| ). They are also in excellent agreement with the 



correlated hyperspherical harmonic (CHH) value of — 24.17±.05 MeV P2| shown by the open 



diamond. The error bar on the CHH value represents the expected truncation error in that 
calculation. However these results are below the Coulomb-corrected Faddeev-Yakubovsky 



value of —24.01 MeV shown by the open circle. The down-pointing triangle shows our 
best VMC upper bound for this case; it is ~ 2 % higher than the exact Eq. 

In the following sections we will give the calculated values of the average energy {Eav) for 
r = 0.04, 0.05, and 0.06 MeV~^ for various A = 6 and 7 states. These provide upper bounds 
to the eigenenergies of the nuclear Hamiltonian H used in this work. The studies of the 
accuracy of GFMC discussed above suggest that, assuming that the low-energy excitations 
in the '^t have been successfully removed by diagonalizing the Hamiltonian matrix in the 
p-shell states, these Eav are at most ~ 0.3 MeV above the eigenenergies for A = 6 states. 
For the ^Li ground state the additional binding obtained by single energy extrapolations is 

47 



only ~ 0.13(15) MeV. This extrapolation was made using a factor five more samples than 
we have for the other states studied, thus no useful extrapolation estimates can be made for 
the other states. We also estimate that the perturbative treatment of if — if' increases Eav 
by less than 0.1 MeV. As will be shown in the next section, the large r behavior of E{t) 
for all the A = 6,7 states is very similar. Thus we expect that the errors estimated for ^Li 
are reasonable approximations to those for other nuclei. By scaling the above two errors 
according to ^^(r = 0) — Eav and {vij) we estimate that the Eav for the A = 7 states is no 
more than ~ 0.5 MeV above the eigenenergies. 

VII. ENERGY RESULTS 
A. Ground States 

The primary results of this paper are the GFMC energies, Eav, of the ten different {J'"]T) 
states in y4 = 6, 7 nuclei shown in Table ^ and in Fig. 0, along with three isobaric analog 
states and the ground states of A = 2 — 4 nuclei. The present results for ^He and ^Li ground 
states and the (3"*';0) excited state in ^Li are all slightly lower, but within error bars, of the 
r-averaged results reported in Ref. The slight improvement may be due to the improved 
while we have obtained much better Monte Carlo statistics than previously, both by 
access to increased computer resources and by more efficient program implementation. The 
r-extrapolated results of Ref. [@ were significantly below the present results, but also had 
a large uncertainty associated with the extrapolation; taking that uncertainty into account, 
the two calculations are consistent. In the particular case of the ^Li ground state, our 
earlier result of —32.4(9) MeV was not inconsistent with the experimental binding of —31.99 
MeV. However, from the various extrapolation tests discussed above, and the propagation to 
r = 0.1 MeV~^, we are now confident that the binding energy with the present Hamiltonian 
is not more than —31.6 MeV. 

Table ^ also gives the VMC energies from the simple starting trial function, and 
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from the more sophisticated \I/y. In A = 3,4 nuclei, picks up about 60% of the energy 
difference between and the final GFMC results. However, for the A = Q,7 nuclei there 
is a much bigger gap between VMC calculations with and the GFMC energies, and the 
\l/y results recover only 20-25% of this energy difference. Clearly there is some important 
aspect to p-shell variational wave functions that is missing from the current ansatz. 

The Argonne f is + Urbana IX Hamiltonian was constructed to reproduce the experi- 
mental binding energies of ^H, '^H, and the equilibrium density of nuclear matter. From 
Table g and Fig. |ll we see that with this Hamiltonian all the A = Q and 7 states studied 



here are underbound. The discrepancy in ^Li and ^Li states is relatively small, < 2% and 
< 5% respectively, and the calculated ground states are stable against breakup into a + d 
and a + 1. On the other hand, the discrepancy in ^He and ''He states is larger, ~ 5% and 
~ 13%, respectively, and the calculated ^He ground state is unstable against a + n + n 
breakup. 

A breakdown of the GFMC energies into kinetic and potential contributions is given 
in Table The kinetic and potential energies grow rapidly as the number of nucleons 



increases, but for a given nucleus, they decrease slightly as the excitation energy increases 
and the nucleus gets more diffuse. The Vijk contribution remains small compared to Vij, 
never exceeding 5%, but because of the large cancellation between K and Vij, it is typically 
25% of the total binding energy. The electromagnetic vjj is dominated by the Coulomb 
interaction between protons, Vci{pp), but about 17% (8%) of its total contribution comes 
from the magnetic moment and other terms in Eqs.( |2.4| - |2.6| ) in He (Li) isotopes. The one- 
pion-exchange term of the potential dominates Vij, providing ~ 70% of the interaction energy, 
while the V^l is smaller than vfj by one order of magnitude. 



Figure |T2| shows the E(t)—E(t=0) from the GFMC calculations of ground states of nuclei 
with 3 < A < 7. The GFMC correction to the VMC (^^t) results has a strong A dependence 
but no significant N — Z dependence. Figure ^ shows [E{r) — E{0)]/\Eav — E{0)\. The 
results for the two s-shell nuclei have the same dependence on r, as do those for the four 
p-shell nuclei. However the p-shell E{t) approach their asymptotic values less rapidly; a fit 
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to the ^Li E{t) for r < 0.03 MeV~^ requires excitation energies of ~ 700 and 90 MeV instead 
of the ~ 1500 and 350 MeV used in the fit for ^He shown in Fig. ^. This is another indication 
that there is a quahtatively new feature necessary for p-shell nuclei which is missing from 
our trial wave functions; this feature does not seem to depend on the N — Z, J'^, or T of 
the nucleus. 



B. Excited states 



A second result of the present paper is the prediction of an additional dozen higher 



excited states obtained in the VMC calculations as shown in Table pCIIj , and in Fig. |14[ We 
have calculated VMC and GFMC excitation energies for eight states, and they agree with 
each other within error bars in all cases. Therefore, we may expect that the VMC excitation 
energies for the other states shown are close to the correct results for this Hamiltonian. Most 
of these higher states are obtained by the diagonalizations within correlated p-shell states 
discussed in Sec. III. 



In Fig. |TT| we see that the difference between the calculated and experimental energies 
increases as A and \N — Z\ increase. However, as seen in Fig. the excitation spectra are 
in good overall agreement with experiment. The states generally occur in the correct order, 
and with reasonable energies. The agreement with the ^Li and low-lying ^Li spectra is very 
good. In particular, the first excited (1^;0) state in ^Li, and the first excited (| ; |) state 
in ^Li, i.e., the first excited states with quantum numbers identical to the ground states, are 
very close to the observed excitations. In the case of ^He and ^He, we predict a number of 
states that have not been observed experimentally, but which could be searched for. A first 
observation of the second (2"^; 1) state in ^He was recently reported the experimenters 



tried to fit their data with single states of different {J'^; T) but did not get a very good fit 
for any one value. Our results suggest there are several states close together in this region, 
which could improve the chances of fitting the data satisfactorily. The states in ^He might 
also be amenable to experimental measurement with the new radioactive beam facilities that 
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are now coming on line. 



C. Isobaric analog states 

Energy differences of isobaric analog states are sensitive probes of the charge- 
independence-breaking parts of the Hamiltonian. To study these it is useful to express 
the energies in an isobaric multiplet, characterized by A and T, in terms of the isospin 
multipole operators of order n: 



EaATz) = E «S:TQn(T,T,) . (7.1) 



V 

n<2T 



The Qn{T,Tz) are orthogonal functions for projecting out isovector, isotensor, and higher- 
order terms [Q; the first terms are Qq = 1, Qi = T^, and Q2 = |(3T| — T^). The coefficients 



^(n) g^j^g then obtained from 



,in) 



aXT = E Qn{T, Tz)EaAT.)/ E Ql^T, T,) . (7.2) 

In first-order perturbation theory, the electromagnetic interaction contributes to the a*^"-* 
for n = 1 and 2, the nuclear CSB potential and kinetic energy contribute to n = 1, and 
the nuclear CD potential contributes to n = 2. The a'-"-' for higher n are zero in first 
order with our Hamiltonian, and there is little experimental evidence for n > ?> terms in 
nuclei |^4[. We have made VMC calculations of the a*-^'^-* in first order by using a CI wave 



function of good isospin, T, and simply varying to compute the Ea,t{Tz)- Table pCIU 



contains results for the T = ^ isovector {n = 1) coefficients in A = 3 and A = 7, the T = 1 
isovector and isotensor (n = 2) coefficients in A = 6, and the T = | isovector and isotensor 
coefficients in A = 7. The energy differences are broken down into v'^, v'"^, v'"^^, and K^^^ 
contributions, with the pp Coulomb, [fci(pp)], other Coulomb (fci?), and magnetic moment 
{vmm) components of v'^ also given. 

The CIB parts of the Hamiltonian induce CD changes in the nuclear wave function, 
leading to higher-order perturbative corrections to the splittings of the isospin multiplets. 
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We have estimated some of these changes in VMC by repeating the calculations using wave 
functions that have a varying Coulomb term, of Eq.( |3.20| ), added to the single-particle 
potential well that is used to generate the (f)p^{Rak) components of \E'y. This results in 
a slightly more diffuse wave function as Z increases, and slightly smaller energy coeffi- 
cients than those obtained with CI wave functions. In the GFMC calculations, the isoscalar 
Coulomb term, f^, provides an additional source of CIB through the propagating Hamilto- 
nian of Eq. (|2.17| ), which depends on T^. However, within the limited propagation time of 



T = 0.06 MeV~^, the main effect of using GFMC wave functions seems to be the slightly 
sharper two-body densities (discussed below) around 1 fm, and consequent changes in the 
CIB potential expectation values. All our GFMC calculations have been made with CD wave 
functions, but a complete set of isobaric analog states was calculated only for A=6, results 



for which are shown on the penultimate line of Table pCIII| . There may also be higher-order 
contributions to the isomultiplet splittings from changes to the CI expectation values, but 
we have no reliable way of extracting these from under the sizeable Monte Carlo errors. 



The results tabulated in Table pClli| indicate that the present Hamiltonian underestimates 



the observed isovector coefficients and gives mixed results for the isotensors. It should be 
remembered, however, that while v*"^ is well determined in the i = partial wave by the 
AW scattering data, it is much less well known in £ = 1 and higher partial waves, while 
the only experimental input for v'"^^ is the nn-pp scattering length difference, which has 
a 20% experimental uncertainty. The A = 3 case would be corrected by a ~ 10% increase 
in v'"^^, but this would not explain much of the discrepancy in the larger nuclei. The 
significant underbinding of the A = 6,7 nuclei with the present Hamiltonian may mean that 
our wave functions for these nuclei are more diffuse than they should be. By far the worst 
discrepancy is for the A, T = 7, | case, where the underbinding of the ground states is also 
the largest. If the Hamiltonian were altered, e.g., by increasing the net attraction from the 
three-nucleon interaction, to obtain the correct binding, the contribution of the CIB forces 
to the isovector coefficients should be increased, both in the dominant pp Coulomb term, 
and the short-range v^^^ . It is more difficult to predict the effect of such changes on the 
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isotensor energy coefficients. 



VIII. ONE- AND TWO-NUCLEON DISTRIBUTIONS 

The one- and two-nucleon distributions of light p-shell nuclei are interesting in a variety of 
experimental settings. For example, the ^He nucleus has been a popular candidate for study 
as a 'halo' nucleus whose last two neutrons are weakly bound. In addition, the polarization 
densities of ^Li and ''Li are important because of possible applications in polarized targets. In 
order to extract information on the spin-dependent nucleon properties from experiments on 
such targets one must, at a minimum, understand the nucleon polarization in the polarized 
nucleus. In this section we provide our results for a variety of nucleon distributions, including 
spin-polarized and averaged single-nucleon densities, spin-dependent and spin-independent 
two-body densities, the proton-proton distributions, and the rms radii, magnetic moments, 
and quadrupole moments. 

As was discussed in Sec. VI, we do not propagate to large enough imaginary time to 
allow the GFMC to significantly modify the rms radius of the p-shell nuclei. Thus they are 
determined almost entirely by the input trial wave function, which is constrained to be near 
the experimental value wherever known. However GFMC does make significant changes to 
densities at small r. A number of the one-body densities are increased near the origin, as 
are the peaks of many of the two-body densities. 



The proton rms radii and static electromagnetic properties are given in Table PCIV . 
These are calculated from using impulse approximation. In general, we know that 
there are significant corrections to the electromagnetic moments from two-body charge and 



current contributions ||45|j46[| . For the magnetic moments, these corrections are only ~ 1-2% 



in isoscalar nuclei like ^H, but are ~ 15-20% in the isovector T = ^ nuclei and ^He. 
Therefore it is not surprising that we see very little discrepancy for the magnetic moment 
in ^Li, but a sizeable error for ^Li; presumably, a calculation including meson-exchange 
contributions would come much closer to the experimental values in the latter case. The 
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quadrupole moment is a more difficult problem, particularly in ^Li, where there is a delicate 
cancellation between the contributions from the deuteron quadrupole moment and the D- 
wave part of the a-d relative wave function. Many cluster models for ^Li fail to obtain the 
observed negative sign; we have trouble getting an accurate measure of the magnitude, for 
reasons discussed above. In the case of ''Li, where there is no such delicate cancellation, the 
value is only ~ 15-20% too low in magnitude. Again, some of this discrepancy might be 
made up by meson-exchange corrections. 

In Fig. |15|, we present the neutron and proton densities for the helium isotopes, calculated 
with GFMC Previously we have found that the ^He charge form factor is in good agreement 



with experimental data in realistic calculations and hence this distribution should 
be quite accurate. As more neutrons are added, the tails of the distributions broaden 
considerably because of the relatively weak binding of the p-shell neutrons. In addition, 
the central neutron and proton densities decrease rather dramatically. This effect does not 
necessarily require any changes to the alpha-particle core, but can be understood at least 
partially from the fact that the alpha particle no longer sits at the center of mass of the 
entire system. The motion relative to the center of mass spreads out the mass distribution 
relative to that of "^He. 

We also find that the small depression obtained in the central density of ^He gradually 
disappears as more nucleons are added. While the depression is clear in the alpha particle, 
it is nearly within our statistical errors in ^He and seems to have disappeared completely 
in ''He. This can again be understood by taking into account the fluctuations in the center 



of mass of the core nucleons about the center of mass of the entire system. Fig. [T^ shows 
the GFMC neutron and proton densities for the lithium isotopes. These densities are very 
smooth functions of the distance from the center of mass. 

The polarization densities for ^Li and ^Li, computed with VMC are presented in Figs. ^ 
and 118], respectively. The spin-up proton density distribution is defined by 
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PvAr) = ^2(^iJ,Mj = J)\Y: l±|ii!ll±|iM5(,_|r,_R_|) |vI;(J,M, = J)) 

(8.1) 

with similar definitions for spin-down protons, spin-up neutrons, etc. Tlie integral of these 
distributions is the total number of spin-up (down) protons (neutrons) in a fully polarized 
state. The integrated quantities can be important in high-energy experiments designed to 
probe the spin- dependence of the neutron or proton structure functions, while their radial 
dependence may be partially accessible in experiments at lower energies. Experimentalists 
are considering using dense, solid polarized ^LiD targets, as an alternative to the deuterated 
ammonia targets currently being used to probe neutron properties. 

Polarization densities have been studied previously in cluster models, with a fixed (unpo- 
larized) alpha core plus interacting valence nucleons. Our calculations include the possibility 
of the spins in the the core alpha particle being polarized by the valence nucleon's spin and 
orbital angular momentum, and hence it is interesting to examine the results for both the 
distributions and the integrated quantities. 

In the spin projection M = 1 state of the deuteron, the polarization of the neutron 
differs from unity because the tensor interaction induces a D-state in the wave function. 
Integrating Pniif) over r yields a probability for up-spin neutrons of 

P(nT) = Ps + \Pd = 1- IPd , (8.2) 

where Ps and Pd are the S- and D-wave probabilities of the deuteron. In the simplest 
two-body (alpha plus deuteron) model of ^Li, the polarization of the neutrons is determined 
by the D-state probability in the deuteron and the D-state probability in the a — d wave 



function. In three-body models, recent calculations [^] have found that the valence neutron 
had P{n^) = 0.93(1). Since in such a model the two core neutrons are unpolarized, this 
corresponds to a total projection P^W]) of 1.93, or a polarization of 29%. 

As expected, the up spins dominate the down spins in the M = 1 state of ^Li at large 
distances from the center of mass. At very large separations, the ratio will be determined 
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solely by the asymptotic D/S state normalization of the alpha-deuteron wave function and 
the D-state probability in the deuteron. At small r, we find that the spin-down density 
exceeds the spin-up density, presumably because the spins of the outer nucleons prefer to 
try to pair with the core nucleons to a spin-zero state. Even though we find this effect to 
be significant, the integrated spin densities agree reasonably well with the cluster model 
calculations. The integrated neutron densities in ^Li are found to be 1.93 for spin up and 
1.07 for spin down, respectively, yielding the same net polarization of 29%. 

The -P(pT) ill ^Li is found to be 1.94 instead of 2 as predicted by the independent- 
particle shell model. The neutrons carry about half this remaining spin, as the spin-up and 
spin-down neutron integrated densities are 1.98 and 2.02. 

We have also computed a variety of two-nucleon distribution functions. These are defined 

by 

P2,p(r) = ^(^1 EOf/(r - |r. - v,m , (8.3) 

i<j 



where the operators O^j are given in Eq.( |3.9| ). While typically these two-body correlations 
cannot be directly measured, they provide the expectation values of two-body operators and 
can be important ingredients in interpreting the results of experiments. Transition matrix 
elements of this type are needed for extracting the effective weak ttNN coupling constant in 



parity-violating experiments. Fig. shows VMC and GFMC calculations of the SijTi ■ tj 
NN pair distribution function. This correlation is strongly dominated by pion exchange, 
and is responsible for the toroidal shapes which characterize the coupling of spin to space 



in the nucleus |^^. We see that the structure is somewhat enhanced by the GFMC. 



Finally, we present the proton-proton distributions (scaled to have normalization inte- 
grals of Z — 1) for ^He, ^He, ^He, ®Li, and ''Li in Fig. ^ These distributions are directly 
related to the Coulomb sum measured in inclusive longitudinal electron scattering; such 
measurements in '^He have been used to put constraints on the ppp{rij), and realistic cal- 
culations agree with the experimental results iSO]. The behavior of Ppp{r) at short distances 



is largely determined by the repulsive core of the NN potential and is nearly independent 

56 



of the nucleus, but at larger distances it is determined by the size of the nucleus. 

We show results for ppp in ^He and ^He in order to directly compare the alpha particle 
proton-proton distribution to that in the alpha-particle cores of ^He and ^He. Unlike the 
one-body densities, these distributions are not sensitive to center of mass effects. We find 
that the proton-proton distribution spreads out slightly with neutron number in the helium 
isotopes, with an increase of the pair rms radius of approximately 4% in going from ''He to 
^He, and 7% to ^He. While this could be interpreted as a sweUing of the alpha core, it might 
also be due to the charge-exchange {ji ■ tj) correlations which can transfer the charge from 
the core to the valence nucleons. Since these correlations are rather long-ranged, they can 
have a significant effect on the proton-proton distribution. 

Finally, we mention that calculations of the electromagnetic ground-state and transition 
form factors are underway. Complete results for these quantities, including exchange current 
effects, will appear in a separate paper. 

IX. CONCLUSIONS 

Quantum Monte Carlo methods are now a powerful tool for the study of light p-shell nu- 
clei. At present, we can write plausible variational wave functions with the proper quantum 
numbers for the given state of interest, but they do not give sufficient binding to provide 
stability against breakup into subclusters. However, the GFMC method rapidly damps out 
the small amount of highly-excited states contained in the VMC wave functions, producing 
ground states that are stable in the case of ^Li and ^Li. The current major hmitation is the 
small value of r that can be reached in most calculations. This makes it important that the 
starting VMC wave functions have a proper mix of the appropriate spatial symmetries, and 
negligible contamination from low-energy excited states. 

The energies obtained for the ground and low-lying excited states are close to, but some- 
what above, the experimental numbers. We believe the discrepancy is probably the fault of 
the Hamiltonian, most likely the phenomenological short-range part of the three-nucleon in- 
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teraction, rather than a shortcoming of the calculation. We note that the difference between 
experimental and theoretical energies is much less that (Vijk), and might be rectified by an 
improved three-nucleon potential. Despite the discrepancies in the ground state energies, the 
excitation spectra are reproduced very well. We believe this is the first demonstration that 
the shell structure of light nuclei can be obtained directly from bare two-nucleon interactions 
that fit AW scattering data. 

The QMC methods developed here can be extended to eight-body nuclei with the present 
generation of computers. We already have calculations in progress for the ground and low- 
lying excited states of A = 8 nuclei. The next major task will be to refine our model 
for the three-nucleon interaction, perhaps including those relativistic corrections which first 
appear at the three-nucleon level [24|,|15|, with the intention of fitting the energies of A=3- 
8 nuclear states with 1% accuracy. Now that accurate QMC calculations of these states 
are possible, there are a host of interesting problems that become accessible, including the 
response of ^Li and ^Li to electron scattering, and many low-energy electroweak capture 
reactions of astrophysical interest, such as ^He((i, 7)^Li and ''Be(p, 7)^B. There also remains 
the problem of adapting the GFMC methods here to the study of larger systems, perhaps 
through methods similar to the cluster-expansion used in VMC calculations of ^'^O E^. 
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TABLES 

TABLE L Values of shape parameters used in generation of trial functions. All units in fm. 
Notation is same as Ref. [17]. 



oo 0.35 Co 

ai 0.4 ci 

at 0.4 d 

ab 0.24 



LI Ro 0.75 

3.0 Ri 2.8 

2.0 Rt 3.7 

Rb 0.4 



TABLE IL Values of asymptotic parameters used in generation of trial functions. Notation is 
same as Ref. [|17[]. 





3H 


A > 4 


^0,0 (MeV) 


3.2 


17.0 


Eo,i (MeV) 


6.0 


16.0 


Eifl (MeV) 


13.0 


23.5 


(MeV) 


6.4 


16.5 


^0,0 


0.85 


1.04 


ho,i 


1.70 


1.71 


hi,o 


1.74 


1.54 


hi,i 


1.72 


1.68 


Vo 


0.026 


0.035 


Vi 


-0.007 


-0.015 


Ci 


0.0003 


0.0003 


at 


0.92 


0.86 
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TABLE III. Values of three-nucleon correlation parameters. Notation is same as Ref. [0. 









-0.00025 








-0.0004 








-0.0005 


y 






0.72 


ql (fm-6) 






0.20 








0.85 








0.16 


qP (fm-l) 






0.05 


qi' (fm-2) 






-0.12 


(72* (fm^^) 






0.12 


9f (fm-2) 






0.85 


ql (fm-i) 






-0.012 


(fm-2) 






0.015 








1.2 


9l 






0.35 


TABLE IV. Energy obtained with different trial functions for 


^Li in MeV, and relative cost to 


compute. 








wave function 


Ey 


dE 


cost 




-25.47(30) 




1.00 




-27.00(13) 


-1.53 


1.19 




-27.11(12) 


-1.64 


1.59 




-27.89(12) 


-2.42 


2.23 




-28.14(11) 


-2.67 


2.66 
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TABLE V. Values of variational parameters in p-shell nuclei. 



^=6 A=7 





i'3S[2] 


1'3d[2] 


1'3P[11] 


2p[3] 


2F[3] 


[21] & [111] 




0.0 


0.0 


0.0 


1.0 


1.0 


1.0 


bsp 


1.0 


1.0 


1.0 


0.0 


0.0 


0.0 


Csp 


0.9 


0.9 


0.9 


0.85 


0.85 


0.85 


dsp (fm) 


3.2 


3.2 


3.2 


3.2 


3.2 


3.2 


dpp 


1.0 


1.0 


1.0 


1.0 


1.0 


1.0 


bpp 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


Cpp 


0.1 


0.1 


0.4 


0.1 


0.1 


0.4 


dpp (fm) 


3.2 


3.2 


3.2 


3.2 


3.2 


3.2 


Rf (fm) 


4.0 


4.0 


4.0 


4.0 


4.0 


4.0 


a/ (fm) 


1.0 


1.0 


1.0 


1.0 


1.0 


1.0 


Vp (MeV) 


-20.0 


-18.0 


-18.0 


-20.0 


-18.0 


-18.0 


Rp (fm) 


4.0 


4.0 


4.0 


4.0 


4.0 


4.0 


dp (fm) 


1.5 


1.5 


1.5 


1.5 


1.5 


1.5 
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TABLE VI. ^LSn components in A = 6 states, listed in order of increasing excitation for T = 

and T = 1. 

(J^;r) iS[2] 3s[2] iD[2] ^j^^] ip[ii] ^^[11] 

(l+,0) — 0.987 — 0.117 0.111 — 

(3+;0) — — — 1 — — 

(2+;0) — — — 1 — — 

(l+,0) — -0.074 — 0.949 -0.306 — 

(l+,0) — -0.153 — 0.300 0.942 — 

(0+;l) 0.967 _ _ _ _ _o.253 

(2+;l) — — 0.880 — — 0.476 

(2+;l) — — -0.476 — — 0.878 
(1+;1) 

(0+;l) 0.262 _ _ _ _ 0.965 



66 



TABLE VII. PlSu components in A = 



7 states, listed in order of increasing excitation for 



r = i and r = |. 



/ TTT 


T) 




F[3J 


2t» roi 1 
^P[21J 








b[lllj 


4S[111] 


(i 


2) 


0.998 




0.001 


0.050 


-0.041 


0.012 






(2 


1 N 
2) 


0.994 




-0.087 


0.001 




-0.068 


-0.010 






1 \ 
2) 




0.998 








0.059 






(i 


2) 




0.995 




0.073 


-0.060 


0.036 






(1 






-0.059 




0.969 


0.168 


-0.171 






(2 


2) 




— U.Uoz 








u.yyy 








2^ 






— n 099 




u.uoy 


U.UIO 






(r 

V2 


2> 


0.035 




0.412 


0.909 




-0.014 


-0.057 




(§" 


2) 






0.864 




0.480 






-0.153 


(r 


2) 






1 












'^2 


2) 










1 








V2 


2) 






-0.448 




0.841 






0.303 
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TABLE VIII. GFMC program performance on IBM SPl and SP2 Wide Nodes. N(A,T) is the 
number of spin-isospin states in the wave function and P is the number of pairs. The columns 
give the times for one propagation step and one energy evaluation and the total time needed for a 
50,000 configuration calculation; see the text for a more complete description. 



N(A,T) P 



Propagation 



Energy Calculation 



msec. 



sec. 



SPl 



SP2 



SPl 



SP2 



Total 
node hours 
SPl SP2 









32 


6 


6.6 


3.4 


0.025 


0.01 


13 


A = 6; 


T = 





320 


15 


63 


26 


0.84 


0.31 


220 


A = 6; 


T = 


1 


576 


15 


100 


40 


1.66 


0.56 


290 


A = 7; 


T = 


1 3 
2' 2 


1792 


21 


460 


170 


10.6 


3.4 


2,230 



7 



TABLE IX. Contributions to the GFMC (0(r)) of Eq. |]8| for ^Li. All quantities are in MeV. 



(0)t 



(0) 



Mixed 



(O) 



Mixed 



{0)t 



(O) 



KCI 



v'c 



^ijk 

Sum 
H' 



143.8(4) 
-168.7(4) 
1.5(0) 
-3.5(1) 
-26.9(1) 
-26.9(1) 



147.3(5) 
475.7(6) 
1.5(0) 
-4.4(1) 
-31.3(1) 
-31.3(1) 



3.5(7) 
-7.0(8) 

0.0(0) 
-0.9(1) 
-4.4(1) 



150.8(10) 
-182.6(11) 
1.5(0) 
-5.4(1) 
-35.7(1) 
-31.3(1) 
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TABLE X. Experimental and quantum Monte Carlo energies oi A = 2 — 7 nuclei in MeV 



'^Z{,r-T) 




VMC (yv) 


KjV mo 


Expt 


2tt/i -I- r»\ 

H(l^; Oj 


-2.2248(5) 






-2.2246 


3x1/1+ 1\ 


-8.15(1) 


-8.32(1) 


-8.47(1) 


-8.48 


4tt f r\-\- r\\ 

*He(0+; 0) 


-26.93(2) 


-27.76(3) 


-28.30(2) 


-28.30 


fijT ^ r\A- 1 \ 

'^He(0^; 1) 


-23.77(6) 


-24.87(7) 


-27.64(14) 


-29.27 


'^He(2^; 1) 


-22.05(6) 


-23.01(7) 


-25.84(11) 


-27.47 


°Li(l^; 0) 


-27.04(3) 


-28.09(7) 


-31.25(11) 


-31.99 


°Li(3^; 0) 


-23.98(7) 


-25.16(7) 


-28.53(32) 


-29.80 


°Li(0^; 1} 


-23.18(6) 


-24.25(7) 


-27.31(15) 


-28.43 


6t ■ /o+ r»\ 

Li(2^; 0) 


-22.58(10) 


-23.86(8) 


-26.82(35) 


-27.68 


fiT> ^ ^ f\A- 1 \ 
r5e(U^; Ij 


-21.73(6) 


-22.79(7) 


-25.52(11) 


-26.92 


He(f ;§) 


-19.02(8) 


-20.43(12) 


-25.16(16) 


-28.82 


7t vS^ 1\ 


-31.59(8) 


-32.78(11) 


-37.44(28) 


-39.24 




-SI 13('8'1 


-32 45n 1 "1 


-36 68 ('30'! 


-38 76 




-25.77(6) 


-27.30(11) 


-31.72(30) 


-34.61 


7Li(f-;l) 


-24.91(7) 


-26.14(11) 


-30.88(35) 


-32.56 


^Li(i";i) 


-18.27(7) 


-19.73(12) 


-24.79(18) 


-28.00 
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TABLE XI. Kinetic and potential energy contributions to GFMC energies in MeV 



^Z{J'^;T) 


K 


Vij 


Vijk 


7 

V- ■ 

V 






2H(1+;0) 


19.81 


-22.05 


0.0 


0.018 


-21.28 


0.0 




50.0(8) 


-57.6(8) 


-1.20(7) 


0.04 


-43.8(2) 


-2.2(1) 


%e(0+;0) 


112.1(8) 


-136.4(8) 


-6.5(1) 


0.86(1) 


-99.4(2) 


-11.8(1) 


6He(0+;l) 


140.3(15) 


-165.9(15) 


-7.2(2) 


0.87(1) 


-109.0(4) 


-13.6(2) 


6He(2+;l) 


131.9(14) 


-155.7(13) 


-7.0(1) 


0.86(1) 


-106.2(5) 


-13.1(2) 


6Li(l+;0) 


150.8(10) 


-180.9(10) 


-7.2(1) 


1.71(1) 


-128.9(5) 


-13.7(3) 


6Li(3+;0) 


146.7(29) 


-174.4(31) 


-7.1(2) 


1.71(2) 


-119.9(5) 


-13.9(4) 


6Li(0+;l) 


135.1(16) 


-161.4(16) 


-6.9(2) 


1.65(1) 


-108.5(4) 


-12.9(2) 


6Li(2+;0) 


139.6(32) 


-166.0(34) 


-6.7(3) 


1.66(3) 


-119.2(5) 


-12.4(4) 


6Be(0+;l) 


134.8(16) 


-160.5(16) 


-6.8(2) 


2.97(2) 


-108.0(4) 


-12.8(2) 


^He(i-;i) 


146.0(17) 


-171.2(17) 


-7.4(2) 


0.86(1) 


-109.9(6) 


-14.1(2) 




186.4(28) 


-222.6(30) 


-8.9(2) 


1.78(2) 


-152.5(7) 


-17.1(4) 




183.0(32) 


-219.1(35) 


-8.2(3) 


1.76(2) 


-151.5(7) 


-16.1(4) 




178.4(28) 


-209.6(30) 


-8.5(3) 


1.78(2) 


-142.2(7) 


-16.1(4) 


7Li(f-;l) 


169.1(31) 


-200.2(33) 


-7.1(3) 


1.73(2) 


-143.2(7) 


-14.2(4) 


^Li(i";i) 


147.8(15) 


-173.8(15) 


-7.2(2) 


1.68(1) 


-109.4(6) 


-13.9(2) 
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TABLE XII. Experimental, VMC and GFMC excitation energies (adjusted to their respective 

ground states) in MeV 

^Z(J^;T) Experiment VMC GFMC 



6He(2+;l) 
6He(2+;l) 
^He(l+;1) 
6He(0+;l) 



1.80 
5.6 

? 
? 



6Li(3+;0) 
6Li(0+;l) 
6Li(2+;0) 
6Li(2+;l) 
6Li(l+;0) 
6Li(l+;0) 



2.19 
3.56 

4.31 
5.37 
5.65 

? 



^Hc(i-;i: 

^Hc(f-;|; 
^He(i-;i: 



^Li(i- 
^Li(|- 
^Li(i- 

^Li(f- 
^Li(|- 

^Li(i- 



I) 

^) 
I) 
\) 
\) 
1) 



0.48 
4.63 
6.68 
7.46 

9.67 

9.90 

? 

11.24 



1.86(10) 
3.61(10) 
3.46(10) 
5.24(11) 



1.80(18) 



2.93(10) 
3.84(10) 
4.23(11) 
5.64(10) 
5.68(11) 
8.96(11) 



2.72(36) 
3.94(23) 
4.43(39) 



0.90(16) 

1.69(16) 
2.08(16) 



0.33(16) 
5.48(16) 
6.64(16) 
9.90(16) 

11.63(16) 
10.14(16) 
10.79(16) 
13.05(16) 



0.76(41) 
5.72(41) 
6.56(45) 



12.65(33) 
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TABLE XIII. Breakdown of VMC isovector and isotensor energy coefficients a^^j, (in MeV) 
obtained with CI wave functions. Total coefficients are given for VMC witfi CI and CD wave 



functions, and for GFMC CD wave functions in ^=6 nuclei. 



A,T,n 


S i 1 


6,1,1 


6,1,2 


7 i 1 

' ' 2 ' 


7 3 1 

' ' 2 ' 


''2' 




0.680(1) 


1.048(2) 


0.186(1) 


1.501(3) 


1.109(4) 


0.119(1) 


[(vciipp))] 


[0.651] 


[1.030] 


[0.167] 


[1.458] 


[1.099] 


[0.114] 


[{VCR...)] 


[0.011] 


[0.014] 


[0.001] 


[0.021] 


[0.012] 


[0.001] 


[(fMM>] 


[0.018] 


[0.004] 


[0.018] 


[0.023] 


[-0.002] 


[0.004] 




0.014 


0.014 





0.025 


0.011 







0.066 


0.035(1) 





0.080(1) 


0.021(2) 













0.101(12) 








0.020(4) 


a^"'y (VMC: CI) 


0.760(1) 


1.097(3) 


0.287(12) 


1.605(4) 


1.141(5) 


0.139(4) 


a^"?^ (VMC: CD) 


0.760(1) 


1.082(3) 


0.277(12) 


1.597(4) 


1.125(5) 


0.132(4) 


a^"!^ (GFMC: CD) 


0.756(1) 


1.120(9) 


0.256(11) 








a^")f, (Expt.) 


0.764 


1.173 


0.223 


1.644 


1.373 


0.175 
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TABLE XIV. VMC values for proton rms radii (in fm), for quadrupole moments (in fm^), and 
magnetic moments (in hn) all in impulse approximation. Only Monte Carlo statistical errors are 
shown; the limitation of GFMC propagation to r = 0.06 MeV— 1 introduces uncertainties of at 
least 0.1 fm in the rms radii. 





VMC experiment 


// 

p" 

VMC 






VMC 


pvTiprnnPTiti 




1.967 


1.953 


0.847 


0.857 


0.270 


0.286 




1.59(1) 


1.60 




2 Q7Q 






^He(i ;i) 


1.74(1) 


1.77 


-1.770(1) 


-2.128 






%e(0+;0) 


1.47(1) 


1.47 










6He(0+;l) 


1.95(1) 












6Li(l+;0) 


2.46(2) 


2.43 


0.828(1) 


0.822 


-0.33(18) 


-0.083 


6Be(0+;l) 


2.96(4) 












^Li(i-;^) 


2.26(1) 


2.27 


2.924(2) 


3.256 


-3.31(29) 


-4.06 


^Be(i-;i) 


2.42(1) 




-1.110(2) 




-5.64(45) 
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FIGURES 

FIG. 1. The experimental spectrum for A = 6 nuclei. 

FIG. 2. The experimental spectrum for A = 7 nuclei. 

FIG. 3. Statistical errors (MeV) in GFMC calculations with 50,000 initial configurations as a 
function of imaginary time. 

FIG. 4. ^He GFMC energy as a function of imaginary time. The fits are described in the text. 

FIG. 5. ^Li GFMC energy as a function of imaginary time. The fits are described in the text. 

FIG. 6. ^Li GFMC energy as a function of imaginary time for various truncations of the 
noncentral parts of ^t- 

FIG. 7. ^He GFMC energy as a function of imaginary time for '^t with various one-body 

FIG. 8. ^Li rms radii as a function of imaginary time for GFMC calculations with '^t of varying 
rms radii. 

FIG. 9. ^Li GFMC energies as a function of the GFMC rms radii from calculations with ^'t of 
varying rms radii. 

FIG. 10. ^He GFMC energy as a function of imaginary time for the Argonne vi4 potential 
without Vijk. Also shown are several previous calculations identified in the text. 

FIG. 11. Spectrum for A=2-7 nuclei from experiment, and in GFMC and VMC calculations. 

FIG. 12. E{t) - E{t = 0) for the ground states of A=^-l nuclei. 
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FIG. 13. E{t) - E{0)/\Eav - E{0)\ for the ground states of A=3-7 nuclei. 

FIG. 14. Excitation spectrum for A=6,7 nuclei from experiment, and in GFMC and VMC 
calculations. 

FIG. 15. The neutron and proton densities in '^He, ^He, and ^He. 
FIG. 16. The neutron and proton densities in ^Li and ^Li. 
FIG. 17. The spin-up and spin-down neutron and proton densities in ^Li. 
FIG. 18. The spin-up and spin-down neutron and proton densities in ^Li. 
FIG. 19. The two-nucleon SijTi ■ tj density of ®Li computed from ^'t and by GFMC. 
FIG. 20. The proton-proton densities in ^He, ^He, ^He, ^Li, and ^Li nuclei. 
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